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ABSTRACT 

The formation and evolution of star cluster populations are related to the galactic 
environment. Cluster formation is governed by processes acting on galactic scales, 
and star cluster disruption is driven by the tidal field. In this paper, we present a 
self-consistent model for the formation and evolution of star cluster populations, for 
which we combine an iV-body/SPH galaxy evolution code with semi-analytic models 
for star cluster evolution. The model includes star formation, feedback, stellar evolu- 
tion, and star cluster disruption by two-body relaxation and tidal shocks. The model 
is validated by a comparison to TV-body simulations of dissolving star clusters. We ap- 
ply the model by simulating a suite of 9 isolated disc galaxies and 24 galaxy mergers. 
The evolutionary histories of individual clusters in these simulations are discussed to 
illustrate how the environment of clusters changes in time and space. It is found that 
the variability of the disruption rate with time and space affects the properties of star 
cluster populations. In isolated disc galaxies, the mean age of the clusters increases 
with galactocentric radius. The combined effect of clusters escaping their dense for- 
mation sites ('cluster migration') and the preferential disruption of clusters residing in 
dense environments ('natural selection') implies that the mean disruption rate of the 
population decreases with cluster age. This affects the slope of the cluster age distri- 
bution, which becomes a function of the star formation rate density (star formation 
rate per unit volume). The evolutionary histories of clusters in a galaxy merger vary 
widely and determine which clusters survive the merger. Clusters that escape into the 
stellar halo experience low disruption rates, while clusters orbiting near the starburst 
region of a merger are disrupted on short timescales due to the high gas density. This 
impacts the age distributions and the locations of the surviving clusters at all times 
during a merger. The paper includes a discussion of potential improvements for the 
model and a brief exploration of possible applications. We conclude that accounting 
for the interplay between the formation, disruption, and orbital histories of clusters en- 
ables a more sophisticated interpretation of observed properties of cluster populations, 
thereby extending the role of cluster populations as tracers of galaxy evolution. 

Key words: galaxies: evolution - galaxies: interactions - galaxies: star clusters - 
galaxies: starburst - galaxies: kinematics and dynamics - galaxies: stellar content 



1 INTRODUCTION 

It is one of the central aims in current astrophysics to 
constrain the formation and evolution of galaxies, and 
their assembly through hierarchical merging (e.g. Sanders & 
Mirabel 1996; Kennicutt 1998a; Cole et al. 2000; Conselice 
et al. 2003; Kauffmann et al. 2003; van Dokkum 2005; Mc- 
Connachie et al. 2009; Hopkins et al. 2010). Galaxy mergers 
play a fundamental role in hierarchical cosmology (White 
& Rees 1978), introducing the evolution of the galaxy pop- 



ulation as a prime tool to verify cosmological models (e.g. 
Kauffmann et al. 1993; Somerville & Primack 1999; Bell 
et al. 2005). Since the late 1980s, observational studies have 
uncovered a wealth of stellar clusters in galaxy interactions. 
Because star clusters are easily observed up to distances of 
several tens of megaparsecs, it is often said that star clus- 
ters can be used to probe the formation and evolution of 
galaxies. This would enable the reconstruction of the merger 
histories of their parent galaxies (Schweizer 1987; Ashman 
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& Zepf 1992; Schweizer & Seitzer 1998; Larsen et al. 2001; 
Bastian et al. 2005). 

The differences between populations of young (massive) 
star clusters and globular cluster systems show the impact 
of nearly a Hubble time of evolution (e.g. Elmcgrccn & 
Efremov 1997; Vesperini 1998, 2001; Fall & Zhang 2001; 
Kruijssen & Portegies Zwart 2009), under the assumption 
that globular clusters initially shared most of the proper- 
ties of current young star clusters (e.g. Elmcgreen & Efre- 
mov 1997). These differences suggest that cluster popula- 
tions can indeed be used to trace galaxy evolution, espe- 
cially because their formation and evolution are known to 
bo governed by their galactic environment (Spitzer 1987; 
Ashman & Zepf 1992; Baumgardt & Makino 2003; Earners 
et al. 2005b; Gieles et al. 2006). It is therefore crucial to 
assess how a cluster population is affected by environmental 
effects. 

There have been substantial efforts in theoretical stud- 
ies to describe the formation and evolution of star cluster 
systems. Possible formation sites of star clusters in general 
and globular clusters in particular have been addressed in 
theoretical studies (e.g. Harris & Pudritz 1994; Elmegreen 
& Efremov 1997; Shapiro et al. 2010) and numerical simu- 
lations (Bekki et al. 2002; Li et al. 2004; Bournaud et al. 
2008; Renaud et aJ. 2008). These studies all point to gas- 
rich environments with high pressures and densities as the 
possible formation sites of rich star cluster systems. How- 
ever, they do not reproduce populations of star clusters and 
globular clusters that are presently observed, because they 
focus on cluster formation and either contain only a very 
simplified description for star cluster evolution or none at 
all. As they age, star clusters leave their primordial regions 
and dynamically decouple from the gas of their formation 
sites. More importantly, star clusters experience extensive 
dynamical evolution after their formation, which shapes the 
characteristics of the star cluster populations that are ob- 
served today. 

Theoretical and numerical studies on the evolution of 
star clusters have shown that clusters dissolve due to two- 
body relaxation in a steady tidal field (e.g. Spitzer 1987; 
Fukushige & Heggic 2000; Portegies Zwart et al. 2001; 
Baumgardt & Makino 2003) and due to heating by tidal 
shocks (e.g. Spitzer 1958; Ostriker et al. 1972; Chernoff et al. 
1986; Spitzer 1987; Aguilar et al. 1988; Chernoff & Wein- 
berg 1990; Kundic & Ostriker 1995; Gnedin & Ostriker 1997; 
Gielcs ct al. 2006). This dynamical evolution leaves a pro- 
nounced imprint on the population that survives disruption. 
In particular, the age and mass distributions of star clus- 
ter populations have emerged as excellent tools to trace the 
disruption histories of clusters (e.g. Vesperini 2001; Fall & 
Zhang 2001; Lamers et al. 2005a; Pricto & Gnedin 2008). 
This implies that the strength of the disruption processes 
will determine how and to what extent the characteristics of 
evolved cluster populations still trace the conditions of their 
formation. 

The census of the formation and evolution of star clus- 
ters has been applied to populations of star clusters in sev- 
eral studies that focus on the modeling of the observed 
cluster age and mass (or luminosity) distributions (e.g. 
Elmegreen & Efremov 1997; Boutloukos & Lamers 2003; 
Hunter et al. 2003; Gieles et al. 2005; Lamers et al. 2005a; 
Smith et al. 2007; Larsen 2009). These studies show that 



the disruption rate of star clusters varies among different 
galaxies, and also that peaks in the age distributions of star 
clusters can be used to trace interaction-induced starbursts. 
Interestingly, the galaxy sample for which the typical disrup- 
tion rates have been derived shows higher disruption rate for 
galaxies with high star formation rates. 

The above analyses of the formation and disruption of 
cluster populations are based on two key assumptions: 

(1) The formation rate of stars and clusters is assumed to 
be constant throughout a galaxy and often also in time. If 
not assumed constant in time, it is paramctcrised with a 
simple function, often a sequence of step functions. 

(2) The disruption rate of star clusters is assumed to be 
constant throughout a galaxy and in time. 

While these assumptions can be made for a first-order ap- 
proach to the formation and evolution of star cluster popu- 
lations, it is known from theoretical principles of star forma- 
tion and cluster disruption that they do not hold in actual 
galaxies. Particularly the observation that the disruption 
rate of star clusters varies among different galaxies shows 
that it should also vary within a galaxy: for individual clus- 
ters as they pass through different galactic regions, but also 
for the entire cluster population as a galaxy evolves. The 
variation with time and space of the cluster formation and 
disruption rates may or may not affect the observable prop- 
erties of star cluster populations. 

Galaxy interactions provide a clear example of the for- 
mation and destruction of cluster populations and the de- 
pendence thereof on the local galactic environment. Large 
numbers of star clusters are formed in the nuclear starbursts 
occurring during galaxy interactions (Holtzman et al. 1992; 
Whitmore & Schweizer 1995; Schweizer et al. 1996; Whit- 
more et al. 1999). These starbursts are triggered by the an- 
gular momentum loss of the gas in the progenitor galaxy 
discs and the subsequent inflow of the gas (Hernquist 1989; 
Mihos & Hernquist 1996; Hopkins et al. 2006). As a result, 
the gas density in certain locations within galaxy mergers 
is very high. It was shown by Gieles et al. (2006) that star 
clusters are efficiently disrupted by the tidal shocks that 
arise when they gravitationally interact with passing giant 
molecular clouds (GMCs). Because the GMC density in cen- 
tral regions of galaxy mergers is high, we should expect tidal 
disruption to be important. This leaves us with the possi- 
ble, intriguing combination of the enhanced formation and 
destruction of star clusters during certain episodes in the hi- 
erarchical buildup of galEixies (see Kruijssen et al. 2010). The 
violent circumstances in the nuclear region contrast with the 
outer regions of a galaxy merger, where cluster formation 
and destruction proceed more temperately. 

In general, the galactic conditions under which star clus- 
ter populations have been formed and have evolved over the 
history of the universe may have varied widely (e.g. Reddy 
& Steidel 2009) . The use of star cluster populations to trace 
galactic histories would therefore benefit from a thorough 
understanding of the co- formation and co-evolution of galax- 
ies and star clusters. Such a census can only be achieved by 
simultaneously assessing all relevant physical mechanisms, 
i.e. combining the formation and evolution of star clusters 
in a single model and allowing for variations with environ- 
ment. 

A thorough way of modeling cluster evolution would be 
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to perform coUisional A^-body simulations of the evolving 
cluster population in a changing galactic tidal field. This 
has been done in several studies, for a limited range of 
cluster masses, orbits and tidal histories (e.g. Vcsperini & 
Heggie 1997; Baumgardt & Makino 2003; Praagman et al. 
2010). However, these papers only consider the effects of 
smooth potentials and ignore the disruptive eflfect of GMCs 
and other substructure in the gas distribution, thereby lim- 
iting the application of such models to globular cluster sys- 
tems and extremely gas-poor galaxies. Moreover, A'^-body 
modeling is computationally so expensive that it is not pos- 
sible to calculate the evolution of the millions of clusters 
that arc formed during the lifetime of a galaxy. In order 
to self-consistcntly compute the formation and evolution of 
an entire cluster population, the best approach would be to 
implement semi-analytic cluster models in numerical simu- 
lations of galaxies. iV-body simulations of dissolving clusters 
in time-dependent tidal fields can then be used to provide 
benchmarks for the semi-analytic cluster evolution models. 

In this paper, we aim to self-consistently model the 
formation and evolution of cluster populations in numeri- 
cal simulations of (interacting) galaxies. For that purpose, 
we have combined semi-analytic star cluster models (SPACE, 
Kruijssen & Lamers 2008; Kruijssen 2009) with an iV- 
body/Smoothed Particle Hydrodynamics (SPH) code for 
modeling galaxies (Stars, Pelupessy et al. 2004; Pelupessy 
2005). As discussed above, the physical mechanisms that 
play a role in the formation and evolution of star clusters 
have all been studied separately in detail. Combining them 
should allow us to obtain a better picture of how different 
galactic environments affect their star cluster populations. 

Pricto & Gnedin (2008) combined cosmological simu- 
lations with a semi-analytic description for the dynamical 
evolution of globular clusters. Their aim was to model the 
population of globular clusters from high redshift to the 
present day, and not the self-consistent modeling of the the 
entire star cluster population including a range of destruc- 
tion and formation mechanisms in different galactic environ- 
ments. Particular examples of differences with our approach 
are the following: 

(1) We include clusters with masses down to a fiducial lower 
mass limit of 100 Mq. In Prieto & Gnedin (2008), only clus- 
ters with initial masses Mi > 10® Mq are considered. While 
this does not influence the present day globular cluster sys- 
tem due to the destruction of lower mass clusters over nearly 
a Hubble time of evolution, it does obstruct the modeling 
of young (globular) cluster populations, in which the cluster 
masses do go down to a physical lower mass limit of around 
100 MQ(see e.g. Portegies Zwart et al. 2010). 

(2) In our simulations, the clusters arc continuously formed 
at the sites of star formation that axe calculated in the N- 
body/SPH simulation. This is one of the main aspects of 
our model that enables a self-consistent modeling of the for- 
mation and evolution of the cluster population. Conversely, 
Prieto & Gnedin (2008) assume that the initial distribution 
of clusters follows the distribution of the baryonic surface 
density, which does not necessarily match sites of star and 
cluster formation. 

(3) In our cluster evolution model, the dynamical mass loss 
rate of clusters due to two-body relaxation depends on en- 
vironmental effects, because it is known that the tidal field 



strength determines the mass loss rate (see e.g. Baumgardt 
& Makino 2003). This aspect of cluster disruption was not 
included by Prieto & Gnedin (2008). 

Smaller differences include the exact prescription for mass 
loss from clusters due to tidal shocks and the evolution of 
the cluster half-mass radius. 

The paper is structured as follows. In Sect. 2, we first 
discuss the simulation code, divided in a section on the 
galaxy (A/^-body/SPH) model, and a section on the deriva- 
tion, construction and testing of the star cluster evolution 
model. The model runs arc summarised in Sect. 3, while 
some key results are presented in Sects. 4 (isolated disc 
galaxies) and 5 (galaxy mergers). The paper is concluded 
with a summary and a discussion of possible improvements 
and potential applications of the model. 

Throughout the paper, we adopt a standard ACDM cos- 
mology and follow the consensus values after the WMAP 
results (e.g. Spergel et al. 2007), with vacuum energy and 
matter densities Qa = 0.7 and Q,m = 0.3, and present-day 
Hubble constant Ho = 70 km s"^ Mpc"^. 



2 SIMULATION CODE 

Our simulations are performed using a coUisionless N- 
body/SPH code in which the formation and evolution of 
star clusters are included as a sub-grid model. The simu- 
lated galaxies contain particles for the stellar, gaseous and 
dark matter components. 



2.1 Galaxy model 

The evolution of the stellar and dark matter components 
arc governed by pure coUisionless Newtonian dynamics, cal- 
culated using the Barnes-Hut tree method (Barnes & Hut 
1986). The particles sample the underlying phase space dis- 
tribution of positions and velocities and are smoothed on 
length scales of approximately 0.2 kpc to maintain the col- 
lisionless dynamics and to reduce the noise in the tidal field 
(which is used for the cluster evolution, see 2.2.2). The Euler 
equations for the gas dynamics are solved using smoothed 
particle hydrodynamics, a Galilean invariant Langrangian 
method for hydrodynamics based on a particle representa- 
tion of the fluid (Monaghan 1992), in the conservative formu- 
lation of Springel & Hernquist (2002). This is supplemented 
with a model for the thermodynamic evolution of the gas 
in order to represent the physics of the interstellar medium 
(ISM). Photo-electric heating of UV radiation from young 
stars is included (assuming optically thin transport of non- 
ionizing photons). The UV field is calculated from stellar 
UV luminosities derived from stellar population synthesis 
models (Bruzual & Chariot 2003). Line cooling from eight 
elements (the main constituents of the ISM H and He as 
well as the elements C, N, O, Ne, Si and Fe) is included. We 
calculate ionization equilibrium including cosmic ray ion- 
ization. Further details of the ISM model can be found in 
Pelupessy et al. (2004) and Pelupessy (2005). 

Star formation is implemented by using a gravitational 
instability criterion based on the local Jeans mass Mj: 

2x3/2 

Mj = ^{^] < M,ef, (1) 
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where p is the local density, s the local sound speed, G the 
gravitational constant and Mref a reference mass scale (cho- 
sen to correspond to observed giant molecular clouds) . This 
selects cold, dense regions for star formation, which then 
form stars on a timcscale Tsf that is set to scale with the 
local free fall time tff. 

= = 

where the delay factor fsi ~ 10. Numerically, the code 
stochastically spawns stellar particles from gas particles 
that are unstable according to Eq. 1 with a probability 
1 — exp (— dt/rsf). The code also assigns a formation time for 
use by the stellar evolution library, and sets the initial stellar 
and cluster population mass distributions (sec below). Me- 
chanical heating of the interstellar medium by stellar winds 
from young stars and supernovae is implemented by means 
of pressure particles (Pclupcssy ct al. 2004; Polupcssy 2005), 
which ensures the strength of feedback is insensitive to nu- 
merical resolution effects. In this way, the global efficiency of 
star formation is determined by the number of young stars 
needed to quench star formation by UV and supernova heat- 
ing, which is set by the cooling properties of the gas and the 
energy input from the stars. 

2.2 Star cluster model 

2. 2. 1 Cluster formation 

Star clusters are formed in the simulations when a Jeans- 
unstable gas particle produces a star particle as described 
above. It is presently not possible to simulate clusters as 
individual particles for the full cluster mass range, because 
even with modern computational facilities it would require 
following too many particles to allow reasonable computa- 
tion times. Therefore, we choose to spawn the star clusters 
as the "sub-grid" content of a new-born star particle. Their 
masses are generated from a power law distribution with an 
exponential truncation at high masses (Schechter 1976): 

NAM oc M"^ exp (-Af/Af,)dA/f, (3) 

where N is the number of clusters, M is the cluster mass, 
and M* = 2.5 x 10® is the adopted exponential trunca- 
tion mass, which is needed to explain the present day mass 
function of Galactic globular clusters (see e.g. Kruijssen & 
Portegies Zwart 2009). We allocate 90% of the new-born 
particle mass for the star clusters (the "cluster formation 
efficiency" or CFE) . Because we adopt a single value of the 
CFE for all particles, its precise value is not important and 
merely acts as a normalisation of the number of clusters. The 
remaining 10% of the mass is considered to be born in un- 
bound associations of stars that immediately disperse into 
the field after they are formed This dispersion does not 
occur physically in the simulation, because the star particle 
contains both the field stars and the star clusters. All stars 
have masses distributed according to a Kroupa (2001) IMF 



^ We do not distinguish between the formation of unbound field 
stars and the early disruption of star clusters during gas expul- 
sion, which is known as "infant mortality" (Lada & Lada 2003; 
Goodwin & Bastian 2006). 



in the mass range 0.08 MQ-rrimax, where rrimax is the maxi- 
mum stellar mass at log (t/yr) = 6.6 (which is the mininmm 
age of the adopted stellar evolution models, see Sect. 2.2.2). 

Owing to the sub-grid nature of the cluster implemen- 
tation, the maximum mass of the formed star clusters is 
determined by the gas particle mass, star formation effi- 
ciency^ and CFE, as the mass of a single cluster can not 
exceed the mass of the star particle it is part of. As a result, 
the typical maximum cluster mass is log (M/Mq) « 5.8. We 
have chosen the number of particles in the sirrmlation such 
as to optimally cover the cluster mass range of interest and 
to have sufficient numerical resolution. An algorithm that 
allows for simultaneous star formation in groups of gas par- 
ticles is being included in a future study. 



2.2.2 Cluster evolution 

After their formation, star clusters evolve by losing mass by 
stellar evolution and dynamical evolution. The star cluster 
evolution is computed with the SPACE models (Kruijssen & 
Lamers 2008; Kruijssen 2009), which include a semi-analytic 
description of the evolution of cluster mass and its stel- 
lar content. They include stellar evolution from the Padova 
models (Marigo et al. 2008), stellar remnant production, 
remnant kick velocities (e.g. Lyne & Lorimer 1994), dynami- 
cal disruption (Lamers et al. 2005a) and the evolution of the 
stellar mass function owing to the stellar mass dependence of 
the ejection rate of stars from the cluster (Kruijssen 2009). 
The total mass loss rate is constituted by the contribution 
from stellar evolution, (dM/dt)ev, and the contribution from 
tidal disruption, [AM / dt)d{a'. 

V dt ; ~ V dt Jev"^ V dt Jdi/ ^ ' 

The mass loss due to stellar evolution is computed using the 
Padova isochrones by following the decrease of the maximum 
stellar mass over one timestep, and integrating the mass 
function within the cluster over the corresponding mass in- 
terval. This method assumes the instantaneous removal of 
stars and ignores the gradual nature of stellar winds. Stars 
typically only lose mass during the last ~ 10% of their life- 
time, during which the maximum stellar mass decreases by 
merely a few percent, and the total cluster mass by even 
less. The mass loss rates of the most massive stars are also 
comparable during the enclosed time interval, which implies 
that the instantaneous removal of the most massive stars is 
balanced by the delay of mass loss from slightly less massive 
stars. This ensures that the obtained mass loss rate is very 
similar to the actual rate due to stellar winds and super- 
novae at any time (see e.g. Kruijssen & Lamers 2008). We 
include the production and retention of stellar remnants as 
discussed in Kruijssen (2009). 

The mass loss rate due to disruption is driven by two- 
body relaxation and tidal shocks: 

/dM\ _ /dM\ /dM\ 

V dt Jdi, ~ V dt Jrix"^ V dt Ah' ^' 



^ This is the fra<;tion of the gas mass that is used to form the 
star particle. 
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where the subscripts 'dis', 'rix', and 'sh' denote disruption, 
two-body relaxation and tidal shocks, respectively. Wc now 
describe the contributions from both mass loss mechanisms^ . 

The mass loss rate due to two-body relaxation is deter- 
mined by the strength of the tidal field and the cluster mass 
(Baumgardt & Makino 2003; Gieles & Baumgardt 2008). 
To describe the mass loss, we adopt the semi-analytic ap- 
proach of Lamers et al. (2005a) that has been extensively 
tested against A'^-body simulations of dissolving star clus- 
ters and observations (e.g. Lamers ct al. 2005b; Giclcs ct al. 
2005; Bastian et al. 2005; Lamers & Gieles 2006). The mass 
decreases exponentially on a disruption timescale that de- 
creases as the cluster mass decreases t^is = (dlnM/dt)~^: 



/dM\ 

V dt /rlx 



M 



(6) 



where the disruption timescale is given by tdis = toM"'. Here, 
the exponent 7 = 0.6 — 0.8 is the mass dependence of the 
disruption timescale, and increases with the King parame- 
ter Wo of the cluster density profile (Baumgardt & Makino 
2003; Lamers et al. 2010). The normalisation constant to is 
the dissolution timescale parameter, which sets the rapidity 
of the disruption and is determined by the tidal field'*. For 
clusters on circular orbits in a logarithmic potential, to has 
been related to the angular frequency of the orbit, and sub- 
sequently to the ambient density Pamb as to oc p~^l^ (Baum- 
gardt & Makino 2003; Lamers et al. 2005b). The physical 
driving force behind cluster disruption is the tidal field. Ac- 
cording to Poisson's law, the tidal field strength T is propor- 
tional to the ambient density, implying a more fundamental 
relation: 



to=to,©(T/T0) 



-1/2 



(7) 



where to,© is the dissolution timescale in a logarithmic po- 
tential at the galactocentric radius of the sun -Rgc,© and 
Tq ^ 7.01 X 10^ Gyr"^ is the tidal field strength at that 
location for a circular velocity of 220 km s~*. For 7 = 0.62 
one obtains to,© = 21.3 Myr, while for 7 = 0.70 we have 
to,© = 10.7 Myr (Kruijssen & Mieske 2009). We adopt a 
density profile with King parameter Wo = 5 for the clus- 
ters and consequently 7 = 0.62. This 'typical' King pa- 
rameter is consistent with observations of open clusters 
(Portegies Zwart et al. 2010) and rapidly dissolving glob- 
ular clusters (McLaughlin & van der Marel 2005; Kruijs- 
sen & Mieske 2009). Clusters with lower King parameters 
{Wo ~ 3) arc susceptible to rapid disruption due to stel- 
lar evolution-induced mass loss (Fukushige & Heggie 1995; 



^ We neglect a third mass loss mechanism, namely the dynamical 
mass loss that is induced by the shrinking of the Jacobi radius 
resulting from the mass loss due to stellar evolution (Lamers et al. 
2010). This is allowed if clusters initially do not fill their Roche 
lobes. In the irregular tidal fields that we axe considering, the Ja- 
cobi radius constantly changes. This implies that the equilibrium 
situation of a cluster filling its Roche lobe is unlikely to occur. 
■* Throughout the paper, wc do not only use the term 'disruption 
time(scale)', but also the more intuitive 'disruption rate', which 
is related to the inverse of the timescale. While the disruption 
timescale is specific to the properties of a cluster and depends on 
its mass and (for tidal shocks) half-mass radius, the term 'dis- 
ruption rate' is used to refer to the general 'disruptiveness' of the 
environment. 



Baumgardt & Makino 2003; Lamers et al. 2010), while high 
King parameters of Wo > 7 are typically achieved after core 
collapse of the most massive systems such as old globular 
clusters. To illustrate the influence of the concentration on 
cluster disruption, we also consider the case of Wo = 7 in 
the rest of the derivation of the model. 

To determine the tidal field strength, we first evaluate 
the tidal field tensor 

where <!> is the gravitational potential and Xi is the i-th com- 
ponent of the position vector. In the simulations, the tidal 
tensor is computed by numerical differentiation of the force 
field, which is smoothed on scales of 200 pc, thereby min- 
imising the sensitivity of the evolution of the star clusters to 
discreteness noise. We use 1% of the smoothing length for 
the differentiation interval. The tidal tensor has three eigen- 
vectors, which denote the principal axes of the action of the 
tidal field. The corresponding eigenvalues represent the mag- 
nitude of the force gradient along these axes, with negative 
eigenvalues denoting comjircssive components of the tidal 
field, and positive eigenvalues indicating extensive compo- 
nents (e.g. Renaud et al. 2008). The tidal field strength T, 
i.e. the quantity that sets the tidal boundary of the cluster, 
is thus equal to the largest eigenvalue of the tidal tensor. If 
the tidal field is fully compressive, i.e. all eigenvalues of the 
tidal tensor are negative, we assume (dM/dt)rix = 0. The 
eigenvalues are computed with the routines by Kopp (2008). 

Tidal shocks disrupt star clusters by increasing the en- 
ergy of the stars that are bound to the cluster. It was shown 
by Kundic & Ostriker (1995) that the first- and second-order 
energy inputs induced by tidal shocks contribute more or less 
equally to the disruption of star clusters, while higher-order 
terms can be neglected. For the mass loss rate due to tidal 
shocks we write 
/dM\ __M 

V dt / sh tsh ' 

where tsh denotes the disruption time for tidal shocks. It 
can be separated in the disruption times due to the first- 
and second-order energy input, tsh,i and tsh, 2: 



(9) 



tsh — (ish,l + *sh,2) 



(10) 

Both components of the disruption time depend on sev- 
eral properties of the cluster and its environment, and will 
change over time. 

The derivation of tsh has been treated extensively in lit- 
erature (e.g. Spitzer 1958, 1987; Ostriker et al. 1972; Kundic 
& Ostriker 1995; Gnedin & Ostriker 1997; Gieles et al. 
2007a; Prieto & Gnedin 2008), though the details of these 
approaches differ. For example, some studies correctly ob- 
serve that not all of the energy input by the tidal shock is 
converted into mass loss (Gieles et al. 2007a), while others 
add the second-order disruption component tsh, 2 (Kundic & 
Ostriker 1995). Also, most studies consider the tidal pertur- 
bation of clusters on closed orbits, for which the frequency 
of disc and bulge shocks is predictable. However, for more 
erratic tidal shocks, tsh should be linked to the tidal field 
(Prieto & Gnedin 2008). Especially when modeling the evo- 
lution of star clusters in galaxy mergers this is an important 
improvement. We follow the lines of most of the above stud- 
ies, and combine their refinements. 
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We first compute the first-order disruption timescale 
due to tidal shocks, which can be expressed as (e.g. Gieles 
et al. 2007a): 



^sh,l — 



At 
f 



E 



AE 



(11) 



where E denotes the cluster energy® per unit mass and At 

is the time since the previous shock. The parameter / is 
the fraction of the relative energy change that is converted 
to a change in cluster mass. This number is smaller than 
unity, because stars escape the cluster with velocities above 
the escape velocity. It is defined as f = din M /din E, and 
has been found to be / ~ 0.25 for two-dimensional (2D) 
shocks^ (Gieles et al. 2006). The internal energy per unit 
cluster mass E is given by 



r]GM 

' 2rh ' 



(12) 



with 7] ~ 0.4 a proportionality constant (e.g. Spitzer 1987), 
G the gravitational constant, and Vh the half- mass radius of 
the cluster. 

We combine the approaches of Gieles et al. (2007a) and 

Prieto & Gnedin (2008) to express the average energy change 
AE of the ensemble of stars in the cluster as a function of 
their average radial position r and the tidal heating param- 
eter Itid- 



1- 



A£=-(Av)2 = -/tidr2. 



(13) 



The tidal heating parameter is written as a function of the 
tidal tensor (Gnedin et al. 1999; Prieto & Gnedin 2008): 

2 

(14) 



in which the integration is performed over the duration of 
the tidal shock for the particular component of the tidal ten- 
sor. The factor ylw.ij (x) represents a parameterised version 
of the Weinberg adiabatic correction (Weinberg 1994a,b,c). 
It is defined for each component of the tidal tensor and de- 
scribes the absorption of the energy injection by the adia- 
batic expansion of the cluster. The adiabatic correction de- 
pends on the product of the average angular frequency of 
the stars within the cluster uj and the timescale Tij of the 
shock for the corresponding component of the tidal tensor 
(Gnedin & Ostriker 1997, 1999): 



= (1 + 



_2 2\-3/2 



(15) 



The value of to is constant when expressed in A'^-body units 
(Heggie & Mathieu 1986; Gieles et al. 2007a), but when 
converted back to physical units it becomes: 



(16) 



^ This is the sum of the internal kinetic energy and the internal 
potential. 

^ Most tidal shocks occur in the orbital plane of the interac- 
tion between the cluster and the perturbing object, e.g. a GMC. 
This corresponds to a 2D shock. A ID shock resembles a head- 
on encounter with the perturbing object, which is relatively rare 
compared to a more distant passage. 



where Cuj denotes a proportionality constant, which for King 
parameters Wo = {5, 7} is Cuj = {0.68, 0.82} (Gieles et al. 
2007a). The timescale of the shock nj is the time interval 
in which the corresponding component of the tidal tensor 
drops by 39%, coinciding with the definition of one standard 
deviation in a Gaussian distribution. 

Substitution of Eqs. 12 — 14 in Eq. 11 now gives the 
disruption timescale due to the first-order effects of tidal 
shocks: 



(17) 



with the ratio r^/r^ — {2,3.5} for King profile parameters 
Wo = {5,7} (Gieles et al. 2006). This equation should be 
complemented with the disruption timescale due the second- 
order energy input, or "shock-induced relaxation" (Kundic 
& Ostriker 1995), which is expressed as 



At 

U.2 = y 



E^ 



{AEf 



(18) 



where (Ai?)^ denotes the stellar ensemble-averaged mean 
square energy change. Following Kundic & Ostriker (1995), 
we write 

1 



{AEY = (vAv)^ = -/tidwV* 



(19) 



The stellar ensemble-average oj'^r^ then follows: 



w2r4 = GM{r)r = CGMrh, 



(20) 

where M(r) represents the mass within radius r, and the 
constant C is defined as 



M{r)r 
Mrh 



(21) 



For King profile parameters Wo = {5, 7} the values are C = 
{0.81,1.03} (M. Gieles, private communication). 

Substitution of Eqs. 12 and 19 in Eq. 18 now gives the 
disruption timescale due to the second-order effects of tidal 
shocks: 



5v'GM , 5v r\ 

tsh,2 - -^3- itid Ai = Y^^U,! 



(22) 



4/C ri 

Using Eq. 10, the disruption time due to the combined first- 
and second order effects of tidal shocks then becomes: 



tah — (ish,l + ^sh,2) 



1 , 12C r-g , , _ ^ , 

J- + — =r Ish,l = Oshlsh.l, 

or] 



(23) 



where for Wo = {5,7} we have C^h = {0.29,0.36}, indicat- 
ing that the contribution of the second-order energy input 
is most important for low-concentration clusters. The mass 
loss due to shocks is applied upon the completion of a shock 
in any of the components of the tidal tensor. Numerically, 
this means that 7tid = unless a shock is completed, i.e. 
one of the components of \Tij\ reaches a minimum that is at 
most 88% of the preceding maximum^. 

We have thus far not defined any prescription for the 
half-mass radius rh- In semi-analytic models that do not 

In a Gaussian distribution, this contrast coincides with the lo- 
cation of 1(7. 
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contain any information regarding the structural evolution 
of star clusters, this is often related to the (initial) mass 
according to a power law relation: 

where rh,4 is the half-mass radius of a 10"' Mq cluster, M(i) 
represents the (initial) cluster mass, and S is the power law 
index. The disruption timescale due to tidal shocks tsh and 
the adiabatic correction Aw both depend on the half-mass 
radius of the cluster (see Eqs. 15 and 17), implying that 
the value of S influences the mass dependence of tsh- It is 
therefore important to include a reliable ijrcscrijjtion for 
the half-mass radius. We have tested several dependences 
of rh on the initial and present cluster mass when compar- 
ing the models to the A'^-body simulations by Baumgardt &: 
Makino (2003). The best agreement is found for S = 0.225 
and rh,4 = 4.35 pc, and when using the present-day mass 
(see Sect. 2.2.3). These parameters are within an accept- 
able range of the 'mass loss-dominated mode' of the radius 
evolution reported in Giclcs ct al. (2011). 

It should be emphasised that the mass-radius relation 
quoted in Eq. 24 does not have the same meaning as the 
mass-radius relation that can be observed for real star clus- 
ter populations (e.g. Larsen 2004). Instead, it approximates 
the evolution of the half-mass radius for a single cluster, 
given a certain initial radius. In a population of star clus- 
ters, which is constituted by clusters of a range of ages, 
initial masses, initial radii, and mass loss histories, the re- 
sulting mass-radius relation of the entire population may 
be very different, as it is set by the collection of states in 
which these different clusters happen to exist at the time of 
the observation. When using a power law formulation, both 
types of mass-radius relation will only be similar if the ini- 
tial half-mass radius of the clusters is also set by the cluster 
mass as in Eq. 24. For mathematical simplicity, we do choose 
to set the initial half-mass radius according to Eq. 24 (see 
Sect. 2.2.3), but it is not a requirement. This approximation 
is supported by clusters in A'^-body simulations, which tend 
towards a well-defined evolutionary sequence (Kiipper et al. 
2008), suggesting that the initial radius may be erased after 
a couple of relaxation times. It is currently not known how 
the half-mass radius evolves in the erratic tidal fields of real 
galaxies, which contain GMCs and spiral arms. We therefore 
choose to adopt the 'conservative' formulation of Eq. 24. In 
a future work, we will include a more sophisticated evolution 
of the half-mass radius. 

The mass loss rates due to two-body relaxation and 
tidal shocks are combined with a model to compute the evo- 
lution of the stellar mass function of the dissolving clusters 
(Kruijssen 2009). In most cases, two-body relaxation gives a 
depletion of the mass function at low masses, because low- 
mass stars have a higher probability to escape than massive 
stars (Hcnon 1969; Vesperini 1997; Takahashi & Portegies 
Zwart 2000; Baumgardt & Makino 2003; Kruijssen 2009). As 
a result, the integrated photometric properties of star clus- 
ters evolve due to dynamical disruption (e.g. Baumgardt &: 
Makino 2003; Kruijssen 2008). By including this, we can use 
the presented models to trace dynamical information of the 
simulated galaxies down to the stellar level. The stars that 



are lost from clusters are added to the field star population 
of the star particle in which the cluster resides. 

The star cluster model is implemented in the galaxy 
evolution code to operate 'on the fly', simultaneously with 
the galactic evolution, rather than having the cluster evo- 
lution calculated a posteriori as in Prieto & Gnedin (2008). 
While this approach is already beneficial because it poten- 
tially allows for a two-way interaction between a galaxy and 
its cluster population, it also implies that the tidal history 
of each cluster is only saved for the most recent time steps, 
which improves the memory efficiency of the simulation and 
allows us to model the full star cluster population all the 
way down to our adopted minimum mass of 100 M0. 

2.2.3 Star duster model testing: method 

We have compared the prescription for star cluster evolu- 
tion from Sect. 2.2.2 to the A^-body models of dissolving 
star clusters by Baumgardt & Makino (2003). These simu- 
lations follow the dynamical evolution of initially Roche-lobe 
filling star clusters in a logarithmic potential with a circular 
velocity of 220 km . The runs contain clusters on circu- 
lar and eccentric orbits between galactocentric radii in the 
range 2.833-15 kpc. The stars in the clusters follow King 
(1966) density profiles with Wo = 5 or Wo = 7, and the 
stellar masses are distributed according to a Kroupa (2001) 
initial mass function between 0.1 and 15 Mq. 

In this section, we exclusively consider clusters on eccen- 
tric orbits, because these clusters experience shocks during 
each pericentre passage. This allows us to test both dis- 
ruption mechanisms rather than only two-body relaxation 
in a steady tidal field, for which the semi-analytic model 
has been tested extensively in previous studies (e.g. Lamers 
et al. 2005a). While the mass loss rate due to two-body re- 
laxation contains no free parameters (see Eqs. 6 and 7), the 
description for tidal shocks depends on the adopted rela- 
tion between the half-mass radius and the cluster mass (see 
Eq. 24), which is governed by the parameters 5 and rh,4- 
Their values are obtained from the comparison. 

To compare our models to the simulations from Baum- 
gardt & Makino (2003), in Fig. 1 we show the time after 
which 95% of the initial cluster mass is lost {tgs%) for our 
models and for their A/^-body runs. The figure shows poor 
agreement if only two-body relaxation is included, but good 
agreement when tidal shocks are accounted for. Addition- 
ally, the infiuence of 6 and rh,4 on tgs^ is shown. As can 
be expected from Eq. 24, 5 affects the mass dependence of 
the disruption time due to shocks (increasing 5 reduces the 
contrast between the disruption times of different masses), 
while rh.4 impacts the normalisation of the disruption time 
(compact clusters live longer). 

As was mentioned in Sect. 2.2.2 and is visible in Fig. 1, 
the best match between our models and the AT-body runs is 
found for S = 0.225 and rh,4 = 4.35 pc. These values should 
therefore approximate the actual evolution of the half-mass 
radii in the AT-body models of the clusters. This is verified in 
Fig. 2, where our adopted mass-radius relation is compared 
to the actual evolution of the half-mass radii of the clusters 
in Fig. 1, showing good agreement. The clusters follow evolu- 
tionary tracks in the mass-radius plane that are very similar 
to each other, indicating that the clusters tend to evolve to 
a common cooling track, analogous to a 'main sequence of 
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Figure 1. Comparison of our analytically estimated cluster lifetimes with the lifetimes found in A'^-body simulations of clusters on 
eccentric orbits (eccentricity e = 0.5) from Baumgardt & Makino (2003). Connected symbols represent different initial cluster masses, 
characterised by different numbers of stars (8k, 16k, 32k, 64k, and 128k), while the solid line shows the 1:1 correspondence between our 
tg^% and that from the A^-body models. Left: without tidal shocks. Middle: Including tidal shocks, for S = {0.15,0.225,0.35} (triangles, 
diamonds, squares) and rji 4 = 4.35 pc. Right: Including tidal shocks, for 5 = 0.225 and rj^ 4 = {3,4.35,5} pc (triangles, diamonds, 
squares). The values at which the agreement between both approaches is best arc written in boldface and are denoted by diamonds in 
the figure. 
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10^ 10* 10^ 10^ 
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Figure 2. Evolution of the half-mass radius as a function of the 
remaining mass for the A'-body simulations shown in Fig. 1, with 
from left to right the initial numbers of stars being 128k, 64k, 32k, 
16k, and 8k (coloured irregular lines). The solid line shows our 
adopted relation between half-mass radius and mass, with rjj 4 = 
4.35 pc and <5 = 0.225. The dashed and dotted lines show the 
variations of rii_4 and 5 from Fig. 1, with dashed denoting rj; 4 = 
{3,5} pc (bottom, top) and dotted denoting 5 = {0.15,0.35} 
(shallow, steep). 

star clusters' as discussed by Ktipper et al. (2008). The ob- 
tained mass-radius relation implies that upon losing stars 
due to disruption, clusters will always slowly evolve towards 
filling their Roche lobes, because the Jacobi radius depends 
on mass as rj oc M^^^ , implying rh/rj oc M~'^'^. The slope 
of the mass-radius relation is also consistent with the 'mass 
loss-dominated mode' from the work by Gieles et al. (2011). 

In principle, the mass-radius evolution of the clusters 
could be a relic of the initial conditions of the A''-body simu- 



lations, in which the clusters initially fill their Roche lobes. 
However, we are considering clusters on eccentric orbits, for 
which the tidal radius continuously changes, suggesting that 
whether or not a cluster initially fills its Roche lobe may be 
irrelevant after a couple of orbits. This would be even more 
important in more realistic, erratic tidal fields. Most impor- 
tantly though, the evolution of the half-mass radius shown 
in Fig. 2 also includes the time after core collapse, when 
any possible imprint of the initial conditions will have been 
erased. Therefore, we do not expect that the details of the 
initial conditions of the A^-body simulations would affect the 
slope or normalisation of the mass-radius relation, especially 
given its simplicity. Nonetheless, the mass-radius relation of 
clusters in erratic tidal fields could deviate from our adopted 
one. We discuss possible improvements of our approach in 
Sect. 6.2. 

2.2.4 Star cluster model testing: numerical resolution 

For any numerical model, it is necessary to check at which 
numerical resolutions the results are reliable. Within a re- 
alistic galactic environment, the tidal field experienced by 
star clusters is very erratic, contrary to the well-defined 
tidal shocks occurring during each pericentre passage in the 
Baumgardt & Makino (2003) simulations. Testing the reso- 
lution requirements of the models (both in time and space) 
should therefore be done for tidal histories that are taken 
from our simulations. 

To explore how the time resolution of the tidal field af- 
fects our modeled star cluster disruption, we computed the 
mass evolution of a 2 x 10'* Mq cluster for 200 tidal histo- 
ries that were randomly drawn from the particles in one of 
our galaxy disc simulations. For each of these histories, the 
evolution was computed seven times, using fixed time steps 
of {1, 2, 4, 8, 16, 32, 64} x 0.932 Myr. For each time step, we 
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time step [yr] 

Figure 3. The mean total disruption time of a 2 x 10"* cluster 

for different time steps, sealed to that for the smallest time step 
(diamonds/solid line). The horizontal dashed line indicates unity, 
which coincides with the leftmost diamond per definition. The 
dark grey tinted area spans the space covered by one standard 
deviation of the distribution at each time step, while the light grey 
represents two standard deviations. The dotted line represents the 
relation for when ttot/ttot, hires were proportional to the time step. 
The maximum time step used in our simulations is indicated by 
the vertical dashed line. 

then scaled the total disruption time of the cluster ttot to 
the total disruption time found for that cluster when using 
the smallest time step (ttot, hires)- This ratio can be used to 
trace the relative change of the total lifetime due to resolu- 
tion effects. Because tidal shocks are events with a certain 
duration, some of them could be skipped when increasing 
the time step, suggesting that in that regime ftot/ttot, hires 
becomes proportional to the time step. 

The mean ttot/itot. hires of the 200 tidal histories is 
shown as a function of the time step in Fig. 3. The relation 
that would be expected if ttot/itot, hires were proportional to 
the time step is also included. The figure shows that for 
large time steps (> 10 Myr), the total disruption time in- 
deed becomes proportional to the time step, as the durations 
of some shocks are then short enough to be skipped, while 
for smaller time steps the total disruption time converges. 
The maximum time step of the particles in our simulations 
(3.73 Myr) is such that time resolution effects should not 
play an important role, particularly because the maximum 
time step is only used for very weaJily accelerated particles 
in dynamically quiet regions. In the simulations, we do not 
use fixed time stops, but adaptive ones instead, increasing 
the resolution as the force on a particle increases (Pclupcssy 
et al. 2004; Pelupessy 2005), up to a maximum resolution 
increase of a factor 4096 (potentially yielding a time step of 
~ 1000 yr). This ensures that tidal shocks, which typically 
occur when the force on a particle is non-negligible, arc al- 
ways well-resolved. In this way, we minimise the effect of the 
time step on our computed cluster lifetimes. 

Whether or not the evolution of star clusters is affected 
by the spatial resolution of the simulations depends on the 
smoothing length and the number of particles used. The 
distribution of mass needs to be resolved in sufficient de- 
tail to ensure that encounters with individual particles do 
not disrupt the clusters. Such disruption would be artificial, 



because individual particles are discrete representations of 
a continuous mass distribution. Whether the resolution re- 
quirements are satisfied can be easily checked with an order- 
of-magnitude estimate. 

Most of the disruption due to an encounter with an 
individual particle would be caused by the correspond- 
ing tidal shock. The presented simulations use a smooth- 
ing length of ft = 200 pc and typical particle masses of 
Mpart = 8 X 10® M0 (see Sect. 3). The typical duration of 
an encounter with a single particle is then approximately 
h/a, with (J the velocity dispersion in a galaxy disc, which 
is of the order 20 km s~^. This gives a typical shock du- 
ration of about 10 Myr. Since we are interested in an up- 
per limit to the disruptive effect of individual particles, we 
ignore the adiabatic correction (Eq. 15) and assume that 
throughout the shock the heating is equal to the tidal heat- 
ing encountered when the cluster is located at the centre 
of the particle. For a spline kernel smoothing, the central 
density of a particle is pcentre ~ Mpart /(ti"/!^). Due to the 
symmetry of a head-on encounter, the tidal tensor is diag- 
onal with values Ttj — — 4GMpart<5ij/(3/i^), which for the 
quoted shock characteristics gives a tidal heating parameter 
of /tid ~ 10^ Gyr~^. If this type of shock would be continu- 
ously repeated over the entire lifetime of a cluster, it would 
take well over 120 Gyr to destroy a 10'* Mq cluster. As is 
evident from Fig. 1 and later sections of this paper, such 
a disruption time is 1-3 orders of magnitude larger than 
typical disruption times. We conclude that for our choice 
of particle numbers and smoothing length, encounters with 
individual particles do not play an important role. Instead, 
the shocks that lead to the disruption of clusters are pro- 
duced by groups of particles, such as spiral arms or com- 
plexes of molecular clouds, which do have a physical mean- 
ing. Consequently, the spatial resolution requirements are 
satisfied. Note that this strongly depends on the smoothing 
length h, because for the maximum tidal heating we have 
7tid oc Tjj- oc Mpart This implies that it is not possi- 
ble to adopt a much smaller smoothing length, which would 
require require vastly larger numbers of particles to reduce 
the particle mass and minimise the effect of encounters with 
individual particles. 

The stability against resolution eflFects is illustrated in 
Fig. 4, which shows the dependence of the star formation 
rate and the number of clusters per unit star formation rate 
on the spatial resolution. The figure shows that the forma- 
tion rates of stars and clusters converge with increasing res- 
olution. The bottom panel gives a measure for star cluster 
disruption, and shows that the variation of the disruption 
rate with spatial resolution is of the order of the inherent 
scatter on the number of clusters. The number of clusters 
very slightly decreases for lower resolutions, because encoun- 
ters with individual particles then become more important 
due to their higher masses. Simulations at higher resolutions 
also exhibit a slight decrease of the number of clusters, be- 
cause in these the structure in the spatial distribution of the 
gas is resolved in more detail. While this may induce a small 
increase of the disruption rate, we do not expect that this 
continues at even higher resolutions, because the amount 
of tidal heating scales with the square of the mass of the 
structure causing the tidal shock, implying that resolving 
increasingly smaller structures results in a correspondingly 
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Figure 4. Dependence of the simulations on the spatial resolu- 
tion. Top: star formation history as a function of time. Bottom: 
number of clusters per unit star formation rate as a function of 
time. Different lines denote simulation IdB (see Sect. 3) run at 
different spatial resolutions. Particle numbers of {1/4,1/2,1,2} 
times those used in simulation IdB are represented by {dash- 
dotted red, dashed green, solid cyan, dotted blue} lines. The 
bumps in the bottom panel for simulations ldBj/2 (^t * = 
3.2 Gyr), IdB (at t = 3.5 Gyr) and ldB2 (at t = 4.1 Gyr) occur 
shortly after the (random) formation of holes in the gas due to 
feedback effects (see text). 



smaller addition to the tidal heating . Figure 4 also illus- 
trates that the mean disruption rate of star clusters is more 
sensitive to random fluctuations than the overall star for- 
mation rate. At certain times, the number of clusters briefly 
increases due to a decrease of the disruption rate. This is 
caused by the random, transient excavation of the gas due to 
feedback in certain star-forming regions, which cause large 
numbers of clusters to experience less disruption. The mean 
disruption rate, which depends on the distribution of the 
gas, then shows larger scatter than the star formation rate, 
which to good approximation is determined by the mean 
surface density of the gas. 



3 SUMMARY OF THE MODEL RUNS 

We construct model disc galaxies with parameters that can 
be related to the outcomes of cosmological ACDM galaxy 
formation models (Mo et al. 1998; Springel et al. 2005). They 
consist of a dark halo with a Hernquist (1990) profile®, an ex- 
ponential stellar disc, a stellar bulge (except for one model) 
and a thin gaseous disc, constructed to be in self gravitating 
equilibrium if evolved autonomously (Springel et al. 2005). 
The disc galaxies are initially set up with 10^-10'' particles 
for the dark matter halo, 22,938-51,250 particles for the stel- 
lar component, and 7,688-25,625 particles for the gas. The 
dark matter haloes have concentration parameters related 
to their total masses and condensation redshifts according 
to Bullock et al. (2001), implying that for a fixed mass the 
halo concentration increases with redshift. The total mass 
is related to the virial velocity V^ir and the Hubble constant 
H{z) at redshift z as Myir = VX/\^^GH{z)]. For aU galax- 
ies, the baryonic disc is constituted by a gaseous and stellar 
component, having a mass fraction rrid = 0.041 of the total 
mass, while the bulge (when included) consists of a stel- 
lar component only, having a mass fraction mb = 0.008 of 
the total mass. These mass fractions are chosen to be con- 
sistent with the fiducial values from recent literature (e.g. 
Springel et al. 2005) and are based on the original constraint 
of 0.03 < md < 0.05 by Mo et al. (1998). The fraction of 
total angular momentum that is constituted by the disc (jd) 
is taken identical to m^. The scale- length of the bulge and 
the vertical scale- length of the disc are 0.2 times the radial 
scale-length of the disc, which is determined by the degree 
of rotation (Mo et al. 1998) through the spin parameter 
A = J\E\/GM'^[^ , in which J is the angular momentum of 
the halo and E its total energy. Table 1 lists the remaining 
properties for the various model runs, i.e. the gas fraction of 
the baryonic disc /gas, the total mass Afvir, the spin param- 
eter A, the number of particles in the different components 
of the model galaxies, and the particle masses of the halo 
particles M^^l° and baryonic particles M^'^H . 

The gas fractions of the galaxy models are chosen to 
cover the range from typical star forming galaxies. Most of 
the total masses represent Milky Way type galaxies, with 
the two exceptions being one half and one tenth of that 
mass, enabling simulations of unequal-mass major and mi- 
nor mergers. The parameter A represents the degree of rota- 
tional support, and is set in accordance with typical spiral 
galaxies in cosmological simulations ((A) = 0.045, Vitvitska 
et al. 2002), except in one case, where we evaluate the in- 
fiuence of the radial scale-length of the disc on the cluster 
population. The number of halo particles is chosen to ensure 
sufficient smoothing of the dark matter halo, and the num- 
ber of stellar disc, stellar bulge and gas particles are chosen 
to minimise the mass difference between the particles and 
alleviate two-body effects. 

The model runs for galaxy mergers are initialised by po- 



° Even for a population of GMCs that follows a power law mass 
distribution with index —2, the total tidal heating would (lin- 
early) increase with GMC mass, despite the many more low-mass 
GMCs than massive ones. 



® This density profile is very similar to profiles found in cosmolog- 
ical simulations (Navarro et al. 1996, 1997). The difference only 
occurs at radii much larger than the scale radius, where the den- 
sity profile of the Hernquist (1990) profile falls of as oc rather 
than r~^. This does not affect the results in this paper, because 
our galaxy merger simulations do not include galaxies on very 
wide orbits. 



Star cluster populations in galaxy simulations 11 



Table 1. Details of the initial conditions for the disc galaxy models. 



ID 


/gas 




z 


A 


.^'halo 


iVgas 


jystar 
disc 


^Btar 
bulge 




part 


Comments 


IdA** 


0.20 


1012 


2 


0.05 


10*^ 


10250 


41000 


10000 




8 X lO"' 


low gas fraction 




0.30 


1012 


2 


0.05 


10*' 


15375 


35875 


10000 


10^ 


8 X lO"' 


standard model 


IdC 


0.50 


1012 


2 


0.05 


10*' 


25625 


25625 


10000 


10'^ 


8 X 10^ 


high gas fraction 


IdD 


0.30 


5x1011 


2 


0.05 


5 X 10"' 


7688 


17938 


5000 


10*5 


8 X 10^ 


half mass 


IdE 


0.30 


1012 


2 


0.05 


10'' 


15375 


35875 





lO'' 


8 X 10^ 


no bulge 


IdF 


0.30 


10" 


2 


0.05 


10*5 


15375 


35875 


10000 


10^ 


8 X 10^ 


low mass 


IdG 


0.30 


1012 


2 


0.10 


10^ 


15375 


35875 


10000 


10^ 


8 X 10^ 


high spin 


IdH 


0.30 


1012 





0.05 


10^ 


15375 


35875 


10000 


10^ 


8 X 10^ 


low concentration 


Idl 


0.30 


1012 


5 


0.05 


10« 


15375 


35875 


10000 


10« 


8 X 10^ 


high concentration 



"In solar masses (Mq). 

''To investigate the relative importance of the two disruption mechanisms, these models are also computed for 
disruption excluding tidal shocks (i.e. only two-body relaxation, 'IdA/B^ix') and for disruption excluding two-body 
relaxation (i.e. only tidal shocks, 'IdA/Bgij')- 

"^This model is also computed for i = {1/4, 1/2, 2} times the number of baryonic particles (i.e. 'IdBj'.) 

'^This model is also computed for a constant disruption parameter to = 2 Myr (see Eq. 6) and no tidal shocks ('IdBflx')- 



Table 2. Details of the initial conditions for the galaxy merger models. 
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sitioning two disc galaxy models on Keplerian parabolic or- 
bital trajectoriesi" with initial separations of approximately 
200 kpc. The actual orbit will decay due to dynamical fric- 
tion, which leads to the merging of the galaxies. The orbital 
geometry of an interaction is characterised by the directions 
of the angular momentum vectors of the two galaxy discs 
and the pericentre distance of the parabolic orbit J?peri. The 
angular momentum vectors of the galaxies are determined 
in spherical coordinates by angles 9 (rotation perpendicular 

1" About 50% of the mergers in cosmological simulations are on 
(near-)paj:abolic orbits (Khochfar & Burkert 2006). 



to the orbital plane) and (rotation in the orbital plane). 
These and other relevant parameters are listed in Table 2, 
where the different model runs are summarised. 

The initial conditions in Table 2 arc divided in three 
categories. The first set of eight runs follow a common con- 
figuration, in which the discs rotate in the orbital plane. 
They are used to test the infiuence of the gas fraction and 
mass ratio of the progenitor discs, and of additional progen- 
itor disc properties such as the presence of a bulge, total 
galaxy mass, and spin parameter (or the radial scale-length 
of the disc). The initial conditions for the six subsequent 
runs are constructed to assess the impact of orbital parame- 
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ters on the star cluster population. We rotate the progenitor 
discs to include retrograde rotation and near-polar orbits, 
which represent the opposite extreme with respect to the 
co-planar configurations of the first eight runs. The effect of 
a wider orbit (a larger pcriccntrc distance) is also considered, 
and a 'random' major merger is also included, in which two 
progenitors with different spin parameters are placed on a 
near-polar, prograde-retrograde orbit. The third group con- 
tains the eight 'random' configurations from Hopkins et al. 
(2009) Barnes 1988), which all have equal probabilities 
of occurring in nature. 

All galaxies are generated without any star clusters, and 
we set t = after 300 Myr of evolution to initialise the 
cluster population. As described in Sects. 2.2.1 and 2.2.2, the 
clusters have masses between 10^ and ~ 10^'* Mq, following 
a Schechter (1976) type initial mass function. The chemical 
composition of the clusters is set to solar metallicity, and we 
assume a King parameter of Wo = 5. 

The properties of the simulated disc galaxies and galaxy 
mergers are not intended to cover all of parameter space, but 
instead should provide a first indication of how the modeled 
star cluster populations are affected by their galactic envi- 
ronment. This sot of simulations represent a basic library 
that can be used to predict certain characteristics of star 
cluster populations and to see how well the simulated clus- 
ter populations compare to observations. 



4 ISOLATED DISC GALAXIES 

As a first application of the model, we consider the simula- 
tions of the isolated disc galaxies from Table 1. As discussed 
in Sect. 2.2.1, the cluster populations are simulated down to 
a lower mass limit of 100 Mq, which for our assumed clus- 
ter formation efficiency and for the typical cluster formation 
and disruption rates of disc galaxies yields about 1 — 3 x 10^ 
clusters per galaxy (also see Kruijssen et al. 2010). Below, we 
discuss the mechanisms driving the evolution of individual 
clusters, and show a number of key properties of the entire 
cluster population. 

4.1 The evolution of individual star clusters in 
disc galaxies 

To illustrate the effects of disruption due to two-body re- 
laxation and tidal shocks, in Fig. 5 we show the orbits and 
the mass evolution of three star clusters with similar ini- 
tial masses {M\ ~ 1.8 x 10* Mq) and times of formation 
{t ~ 2.20 Gyr) from simulation IdB. They arc on differ- 
ent orbits and therefore experience differing tidal evolution. 
Clusters orbiting at small galactocentric radii evolve in a 
stronger tidal field (and generally have smaller Jacobi radii) 
than clusters orbiting at large galactocentric radii. Also, the 
number and intensity of tidal shocks is typically larger for 
clusters orbiting close to the galactic centre, due to the 
higher gas density in their environment. These differences 
result in contrasting mass loss histories and total disruption 
times, as the innermost cluster survives for about 100 Myr, 
the middle cluster persists for about 300 Myr, and the outer 
cluster is disrupted after 400 Myr, even though it has the 
lowest initial mass of the three clusters. The jumps in the 
mass loss history indicate the effect of tidal shocks, which 



are generally stronger (and potentially more frequent) for 
clusters on narrow orbits. Despite these trends. Fig. 5 also 
shows that a cluster can be disrupted by a single tidal shock 
almost anywhere in the galaxy, even at radii well beyond 
the solar galactocentric radius. This was also discussed by 
Gieles et al. (2006), who showed that clusters with masses 
< 2 X 10* Mq can be disrupted during a single encounter 
with a spatially extended CMC of 10® Mq. Relative to dis- 
ruption due to subsequent, small encounters, disruption by a 
single, violent GMC encounter is most prominent for cluster 
masses of about 10* Mq (Gieles et al. 2006). The clusters 
in Fig. 5 are characteristic examples of this. In general, the 
strongest tidal shocks occur at times when the clusters cross 
regions of high gas density, for instance during spiral arm 
passages. This is best seen in the snapshots at t = 2.33 Gyr 
and t = 2 A3 Gyr, because between those snapshots the outer 
cluster is overtaken by a dense region^^, causing it to lose 
almost 75% of its mass due to the rapid change of the tidal 
field. 

Figure 6 illustrates the relation between the mass loss 
and the tidal field for the two long-lived clusters from Fig. 5. 
It shows the evolution of the cluster mass, together with the 
tidal field strength as defined in Sect. 2.2.2, and the running 
integral that is used to compute the tidal heating parameter 
/tid in Eq. 14, which is defined as: 

where tiast is the time of the last shock and t is the cur- 
rent time. This quantity represents the accumulated tidal 
heating since the last shock, and resets after each shock is 
completed. Contrary to Jtid, it docs not include the adiabatic 
correction. Therefore, it is only a measure for the tidal shock 
heating imposed by the tidal field and does not contain any 
information on the response of the cluster experiencing the 
shock. 

The evolution of the tidal field strength for both clusters 
in Fig. 6 shows that it is indeed larger for the cluster orbiting 
at the smaller galactocentric radius. A comparison of the 
tidal field strengths experienced by the clusters just after 
their formation explains why the outer cluster loses its mass 
more slowly initially. The moments at which both clusters 
suffer their first violent mass decrease can be Eissociated with 
jumps in the amount of shock heating, indicating the effect 
of tidal shocks. In the case of the inner cluster, this gives 
rise to its total disruption. The second moment of violent 
mass loss of the outer cluster cannot be coupled with an 
increase of the shock heating, because the shock took place 
in between two snapshots and is therefore skipped by the 
output of the simulation. 

The evolution of the clusters in Figs. 5 and 6 illustrates 
that the disruption rate varies with time for individual clus- 
ters and varies with space when considering the cluster pop- 
ulation as a whole. This is very important when interpreting 
observed star cluster populations. The time- variation of the 
disruption rate for individual clusters can mask the effect of 
a different mass dependence of star cluster disruption. For 



The outer cluster is situated beyond the co-rotation radius of 
the galaxy. 
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Figure 5. Evolution of the orbits and masses of tiiree clusters orbiting at different galactocentric radii in isolated disc galaxy simulation 
IdB. From top to bottom, the consecutive panels show the situation at different times t, while from left to right the orbital evolution 
in the x-y plane (face-on), the orbital evolution in the x-z plane (edge-on), and the mass evolution are shown. The innermost cluster 
is represented by the dark red solid lines, the middle cluster by red dotted lines, and the outermost cluster by blue dashed lines. If at 
any particular snapshot a cluster is still undisrupted, its position and mass are marked with thick dots. The orbital trajectories remain 
visible after the clusters are disrupted. The small grey dots in the x-y and x-z plane views map the distribution of the gas particles in 
the simulation. 
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Figure 6. Evolution of the cluster mass and the tidal field for 
the two outer clusters from Fig. 5, indicated by the same colours 
and line styles. The diamonds in the middle and bottom panel 
indicate the times of each snapshot. Top: The mass evolution. 
Middle: Evolution of the tidal field strength experienced by each 
cluster, defined as the largest eigenvalue of Eq. 8 (see Sect. 2.2.2). 
Bottom: Running integral of the total amount of shock heating 
experienced by the cluster (see text and Eq. 25). 



instance, if the disruption parameter to in Eq. 6 were to in- 
crease with time because a cluster is leaving a region with 
a high gas density (also see Elmegreen & Hunter 2010), for 
that time interval one would derive a lower value of 7, which 
is the mass-dependence of the disruption timescale. 

When considering the space- variation of the disruption 
rate throughout a population of clusters, it is inevitable that 
the mean disruption rate decreases with age^^ , because clus- 
ters in disruptive environments have the shortest lifetimes 
and never reach old ages, while clusters in less violent set- 



Provided that the galaxy and formation sites of the clusters 
do not change much on timescales of ~ 1 Gyr. Galaxy mergers 
are a clear exception to this. 
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Figure 7. Mean age of star clusters (t) as a function of galacto- 
centric radius Rgc for two galaxy simulations at t = 2.5 Gyr. The 
red solid line shows the relation for the physically computed dis- 
ruption rate from Sect. 2.2.2, while the blue dotted line represents 
the result for a simulation with a fixed disruption rate. 



tings become older. This process can be regarded as a form 
of 'natural selection' acting on the star cluster population, 
and tends to flatten the age distributions of star clusters (see 
Sect. 4.4). 

The time-variation of the disruption rate is particu- 
larly interesting in view of recent discussions in literature, 
in which it is debated whether or not star cluster disruption 
depends on the cluster mass (e.g. Fall et al. 2005; Whit- 
more et al. 2007; Gieles & Bastian 2008; Larsen 2009; Bas- 
tian et al. 2009). It is crucial that environmental depen- 
dences are taken into account before inferring any conclu- 
sions about the mechanisms driving cluster disruption from 
observations, because the age and mass distributions of clus- 
ters are susceptible to variations in the environment. This 
holds particular relevance in non-equillibrium settings such 
as interacting galaxies (see Sect. 5). 



4.2 The variation of the disruption rate with 
galactocentric radius 

A second consequence of the variability of the disruption 
rate is related to its variation with space. It implies that the 
properties of the star cluster population, such as the slope 
of the age and mass distribution will depend on the local 
environment within a galaxy. Galaxy- wide distributions may 
indicate the average properties of the cluster population, but 
interpreting them can yield systematic errors when assuming 
the disruption rate is the same everywhere in the galaxy. 
For instance, the effects of cluster disruption are stronger 
towards the galactic centre than in the outskirts of a galaxy, 
implying that the properties of the cluster populations in 
both regions will differ. 

The environmental dependences on the star cluster pop- 
ulation can be qualitatively illustrated by considering two 
variations of simulation IdB. Figure 7 shows the mean clus- 
ter age as a function of galactocentric radius for two galaxy 
disc simulations: one in which the disruption rate is calcu- 
lated as described in Sect. 2.2.2 (model IdB), and one with 
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a disruption rate that is constant in time and space (using 
to = 2 Myr, see Eq. 6) and does not include tidal shocks 
(model IdBflx). For the galaxy with the physically moti- 
vated disruption rate, the spatial distribution of the mean 
cluster age is as expected. The youngest clusters are found 
in the galactic centre, where the star formation rate den- 
sity is highest. Due to the high gas density, clusters in the 
galactic centre disrupt on shorter timescales than in the out- 
skirts of the galaxy, resulting in a mean cluster age that in- 
creases with galactocentric radius. This contrasts with the 
age profile of the cluster population in the galaxy with a 
fixed disruption rate, where the mean cluster age is approx- 
imately constant throughout the galaxy and the scatter is 
strictly due to local variations in the star formation history 
and stochastical effects^^. Observations of cluster popula- 
tions in real galaxies (e.g. van den Bergh & McClure 1980; 
Gieles et al. 2005; Froebrich et al. 2010) show that the mean 
age increases with galactocentric radius, contrary to the re- 
sult for a fixed disruption rate. For the inner disc of M51, 
Gieles et al. (2005) find that the disruption rate varies by a 
factor 1.8 between radial intervals of 1-3 kpc and 3-5 kpc. 
Assuming that the mean age is directly proportional to the 
disruption timescale, this is of the same order of magnitude 
as for the model shown in Fig. 7, for which we find that the 
ratio between the mean ages of the clusters in these intervals 
is 1.4. 

These results substantiate that star cluster disruption 
is indeed driven by environmental effects. Additionally, they 
show that the suggestion of a disruption rate that increases 
with the star formation rate, which is found when consid- 
ering variations between different galaxies (e.g. Boutloukos 
& Lamers 2003; Lamers et al. 2005b), also holds within a 
single galaxy. This is easily understood by noting that both 
the formation and disruption of clusters peak in dense envi- 
ronments. 

4.3 The relative importance of tidal shocks and 
two-body relaxation 

The relative contributions to star cluster disruption of two- 
body relaxation and tidal shocks can be quantified by con- 
sidering the number of clusters in simulations for which ei- 
ther mechanism is neglected. The fraction of the total dis- 
ruption contributed by tidal shocks fsh is then given by 

U^l-^, (26) 

where A^both is the number of clusters in a simulation in- 
cluding both disruption mechanisms (e.g. IdA and IdB), 

The spatial variation of the star formation rate (SFR) cannot 
produce the behaviour of the mean cluster age that is shown in 
Fig. 7. Without a time-variation, the cluster age distributions 
at different galactocentric radii would still yield the same mean 
age, irrespective of the relative formation rates. Within a stable 
disc galaxy, the relative time variations at different galactocentric 
radii are not sufficiently large nor persistent enough to cause a 
spatial trend of the mean cluster age, as is also shown by the line 
denoting the galaxy with a fixed disruption rate. Indeed, Fig. 7 
could have been made at any time in our simulation other than 
the time shown, and the mean age would have shown the same 
spatial variation. 
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Figure 8. Fraction of disruption due to tidal shocks as a function 
of time. The red solid line shows the evolution for all clusters, 
while the light grey dotted line only includes clusters within a 
galactocentric radius of 10 kpc, and the dark grey dashed line 
represents the evolution for the clusters beyond 10 kpc. Top: disc 
galaxy with a gas fraction /gas = 0.20 (simulation IdA). Bottom: 
disc galaxy with a gas fraction /gas = 0.30 (simulation IdB). 

and Arix is the number of clusters in a simulation for which 
only disruption by two-body relaxation is included and tidal 
shocks are not considered (e.g. IdA^ix and IdB^ix). Per def- 
inition Aboth < Arix. For different radial bins, /sh is shown 
as a function of time in the top panel of Fig. 8 for two 
galaxies with different gas fractions (IdA and IdB). The 
contribution by tidal shocks is typically 80-85% of all dis- 
ruption, which is very similar to the analytical estimate by 
Lamers & Gieles (2006) for the solar neighbourhood. The 
value increases with the gas fraction of the disc, because 
GMCs and spiral arms are the most important sources of 
tidal shocks. For relatively gas-rich discs such as in simula- 
tion IdB (/gas ~ 0.30), the contribution from shocks does 
not vary much with galactocentric radius, but for gas-poorer 
discs, shocks are more important in the inner regions of the 
disc. This occurs because beyond a certain galactocentric 
radius the gas density becomes too low to sustain star for- 
mation (e.g Kennicutt 1989; Schaye 2004; Pelupessy et al. 
2004), yielding less or no energy injection by feedback and a 
less filamentary gas distribution, which in turn implies that 
tidal shocks are less important. The characteristic radius for 
this transition is smaller in gas-poor galaxies, which is illus- 
trated by the contrast between the inner and outer parts of 
the disc in the upper panel Fig. 8. Because discs also be- 
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come more gas-poor as they age, the relative contribution 
by shocks slightly decreases with time. 

For the adopted mass-radius relation, the ratio between 
the disruption timescales due to two-body relaxation and 
tidal shocks is trix/ish oc Af^ '^/tid- If this ratio is larger than 
unity, tidal shocks dominate cluster disruption, while a ra- 
tio below unity implies that disruption is mainly driven by 
two-body relaxation. The relative importance of tidal shocks 
increases with cluster mass until a few times 10^ M© , when 
the adiabatic correction in the tidal heating parameter 7tid 
(see Eq. 14) becomes non- negligible and inhibits disruption 
due to tidal shocks. This means that the relative importance 
of tidal shocks peaks at a certain cluster mass. For the pa- 
rameters in this paper, this is /sh ~ 0.9 at M ~ 10* M©, 
but the precise value depends on the mass-radius relation. 

4.4 The age distributions of star clusters in disc 
galaxies 

The balance between cluster formation and destruction gives 
rise to a cluster population with a certain age distribution. 
The age distribution of star clusters is often used as a probe 
to study star cluster disruption (e.g. Gieles et al. 2005; Chan- 
dar et al. 2006) , or to assess the formation history of a galaxy 
(Hunter et al. 2003; Gieles et al. 2005; Smith et al. 2007). 
In order to obtain a reliable interpretation of the cluster age 
distribution, it is important to understand its evolution in 
different galaxies. 

To investigate possible correlations between the cluster 
age distribution and galaxy properties, we have fitted the 
logarithmic slope a of the age distribution {AN /At oc r") 
in the age range log(T/yr) = 7.7-9 for all snapshots of 
our galaxy disc simulations. When constructing the age dis- 
tribution, we consider all available clusters, implying that 
the samples are mass-limited with M ^ 100 Mq. We have 
also fitted the logarithmic slope 13 of the SFR-corrected age 
distributions in that range {[AN /At\/SFR oc t^). The age 
range has been chosen such that the effects of gas expulsion 
due to supernovae are no longer relevant and a sufficiently 
large part of the age distribution is covered to obtain a re- 
liable slope. The fits have been made with 13 bins in the 
specified age range, using a variable bin width to accom- 
modate equal numbers of clusters in each bin. The clusters 
outside the fitted age range are binned using the same num- 
ber of clusters per bin. We have adopted Poissonian errors 
for AN /At, scaling the square root of the number of parent 
star particles instead of using the number of clusters in each 
bin, because the ages of the clusters within a single star par- 
ticle are identical (see Sect. 2.2.1). In practice, this means 
that the relative error decreases with age, because the mean 
number of clusters per particle decreases. To ensure a reli- 
able fit, the slopes have only been measured at times when 
a galaxy contains clusters older than 1 Gyr. 

Given the time interval between subsequent output 
snapshots, the above procedure results in 1175 fitted age 
distribution slopes, covering eight different disc galaxy mod- 
els. These slopes should be considered 'mean' slopes for 
the specified age range, because the age distribution does 
not always follow a single logarithmic slope over the fitted 
age range. This is illustrated in Fig. 9, which shows (SFR- 
(un)corrected) age distributions for two difi^erent galaxies, 
of which the upper one (IdB) is indeed ill- fitted by a sin- 
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Figure 9. Top: Age distributions of clusters in simulations IdB 
(red solid line) and IdG (blue dashed line, vertically offset by 
1.5 dex) at t = 3 Gyr for ages in the range log(T/yr) = 7.7- 
9. The dotted lines represent power law fits to the data in the 
age range indicated by the shaded area. Bottom: Same as above, 
but with the age distributions divided by the star formation rate 
(SFR-corrected). The error bars are computed as described in the 
text. In both panels the fitted slopes are indicated. 

gle power law. For the SFR-corrected age distributions, the 
negative slope is solely the result of disruption, with small 
variations due to stochastical effects. The SFR does not vary 
much in isolated galaxies, and therefore only affects the fit- 
ted slopes by a few hundredths. 

In models with the same disruption rate for all clus- 
ters and a constant SFR, the 'classical' age distribution is 
characterised by two components (e.g. Boutloukos & Lamers 
2003; Lamers et al. 2005a). At young ages, the age distribu- 
tion is flat, because no clusters are disrupted within such a 
short time interval. Beyond the lifetime of the lowest mass 
cluster, the age distribution steepens. This is the disrup- 
tive (old) end of the age distribution, which has a slope of 
P — — 1/7, where 7 is the mass dependence of the disruption 
timescale (see Eq. 6)^*. 

The models presented in this paper assume a more real- 
istic formulation, in which the disruption rate is affected by 
the variation of the tidal held strength and by tidal shocks. 
Nonetheless, for the sake of illustration it is important to 



^ ft is assumed that the logarithmic slope of the cluster initial 
mass function {dN /dAl) is —2. 
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Figure 10. Schematic representation of the two processes leading 
to a cluster disruption rate that decreases with age. Top: Cluster 
migration. Bottom: Natural selection. The large dots mark star 
clusters, the small dots represent debris from disrupted clusters, 
and the clouds denote gas clouds. Time increases from left to right 
in both sequences. 



indicate what the disruption-dominated slope of the age dis- 
tribution would be if the tidal field strength would be the 
same throughout a galaxy, and all clusters would experi- 
ence the same tidal shocks. In such a scenario, the adopted 
value of 7 = 0.62 for disruption due to two-body relaxation 
gives /3 = —1.61. For rapid shocks (i.e. a negligible adiabatic 
expansion during the shock), our adopted mass-radius rela- 
tion yields /3 = —3.08, while for slow shocks (i.e. a dominant 
adiabatic expansion during the shock) we have /? — —1.23. 
Fast disruption (i.e. rapid shocks) thus yields a steeper age 
distribution than slow disruption. 

Evidently, the bulk properties of the cluster population 
are determined by a combination of the above mechanisms, 
covering a range of tidal field and shock strengths. As such, 
the fitted slope of the age distribution is not only determined 
by the mass dependence of the disruption timescale, but also 
by possible trends of the disruption rate with cluster age and 
by the rapidity of disruption in general. 

The difference between the age distributions shown in 
Fig. 9 should be the result of the differences between the ini- 
tial conditions of both simulations. Galaxy IdG only diflers 
from IdB by its (larger) spin parameter (see Table 1), im- 
plying a correspondingly larger scale radius and lower (gas) 
density, which yields a lower disruption rate (see Sects. 4.1 
and 4.2). Because of the more rapid disruption in simulation 
IdB only very few clusters survive for ~ 1 Gyr, causing the 
depletion in the oldest bin, which in turn steepens the fitted 
slope. Again, faster disruption implies less survivors and a 
potentially steeper fitted slope than for slow disruption. 

Another effect is that the disruption rate due to tidal 
shocks will typically decrease as clusters age. This happens 
for two reasons (also see Fig. 10): 

(1) 'Cluster migration': because clusters move out of their 
primordial environment, the ambient gas density typically 
decreases as they age, giving rise to fewer tidal shocks and a 
lower disruption rate at older ages (see Elmegreen & Hunter 
2010). This evolution of the mean disruption rate is more 
pronounced if the density contrast between the star forming 
region and its surroundings is large. 



(2) 'Natural selection': at any given time, clusters in regions 
with a high disruption rate are less likely to survive than 
clusters in low disruption rate regions. Such selection implies 
that at older ages only the clusters in low disruption rate 
regions axe left, causing the disruption rate to decrease with 
age (also see Sect. 4.1 and Fig. 5). This evolution of the 
mean disruption rate is more pronounced if there is a large 
spread in disruption rates, like in galaxies with large density 
contrasts between different regions. 

These two effects make the disruptive end of the age dis- 
tribution shallower and steepen the young end of the age 
distribution. In the extreme case, this can lead to an age 
distribution following a single power law with a slope of — 1 
over the majority of the age range. The effects of cluster mi- 
gration and natural selection are strongest for galaxies with 
low gas densities, because in those galaxies the density con- 
trast between star forming regions and their surroundings 
is larger than in high gas density galaxies per definition^^. 
While already present in simulation IdG, it could thus be 
even more important in dwarf galaxies, which have very low 
gas densities. For out-of-equillibrium systems such as galaxy 
mergers, the dependence of the age distribution on the mean 
gas density is different (see Sect. 5). 

Above, we discussed: (1) the different disruption pro- 
cesses shaping the age distribution, (2) the effect of the 
largest possible cluster lifetime on the fitted slope, (3) the 
effect of cluster migration, and (4) the effect of natural selec- 
tion. For all four of those, galactic environments with high 
gas densities steepen the slope. As discussed at length be- 
fore, cluster disruption is governed by the gas density (pgas), 
implying that in isolated disc galaxies, the fitted slope of 
the age distribution can be used as a measure for the ra- 
pidity of cluster disruption. The gas density also sets the 
star formation rate density (psfr) of a galaxy through the 
Schmidt-Kennicutt law (Schmidt 1959; Kennicutt 1998b). 
One would therefore expect a correlation between the fit- 
ted slopes of the age distributions in different galaxies and 
their star formation rate density. To obtain a measure for 
the star formation rate density that can be determined and 
compared for disc galaxies as well as for galaxy mergers at 
any time during their interaction, we define the mean star 
formation rate density within a sphere with a radius equal 
to the half-mass radius of the gas -Rh.gas: 



Ph,SFR 



SFRh 3 SFRh 



(27) 



with Vh.gas the volume of the sphere, and SFRh the star 
formation rate within Vh,gas- For isolated disc galaxies, most 
if not all of the star formation occurs within i?h,gas- 

We show the relations between Ph.SFR and the fitted 
slope of the cluster age distribution a and fitted slope of 
the SFR-corrected age distribution /3 in Fig. 11 for all 1175 
fits. As expected, it shows an inverse correlation between the 
slope of the cluster age distribution and the star formation 
rate density. For the uncorrected slope a, the fitted relation 
is given by 

Q = C-0.601ogph,SFR, (28) 
where C = —4.66 is a fitting constant that has no particular 

This holds for isolated galaxies. 
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we present a similar result for galaxy mergers, in which the 
number of clusters decreases during a merger despite the 
large starbursts and corresponding cluster production. The 
net destruction of clusters is attributed to enhanced cluster 
disruption that is driven by the high gas density. The analy- 
sis of Kruijssen et al. (2010) does not rely on the cluster age 
distributions, but instead considers the number of clusters 
as a function of time. The number of surviving clusters is 
found to decrease with increasing starburst intensity, which 
is similar to the relation presented here. 

The scatter around the relation between the slope of 
the age distribution and the star formation rate density is 
substantial. Within a single galaxy, a and /3 vary by 0.5 at 
a given star formation rate density, depending on the mo- 
ment at which the galaxy is observed. Because it is relatively 
isolated in the displayed plane, galaxy IdG in Fig. 11 pro- 
vides a clear illustration of the spread. Recent debates in 
literature about the mass-dependence of cluster disruption 
involve differences of a similar magnitude, quoting slopes of 
-1 (Whitmore et al. 2007; Chandar et al. 2010) to -1.5 
(Boutloukos & Lamers 2003; Silva- Villa & Larsen 2010). As 
is shown by Fig. 11, such variations may occur even within 
a single galaxy. Figure 11 also illustrates that a slope of — 1 
is more likely to occur in galaxies with low star formation 
rate densities. As such, both sides of the debate show cases 
that can arise in the framework for star cluster disruption 
that is presented in this paper^^. 



Figure 11. Relation between the fitted logarithmic slope of the 
cluster age distribution in the age range log(T/yr) = 7.7—9 and 
the mean star formation rate density py^ sfr, which is defined for 
a sphere with a radius equal to the half-mass radius of the gas. 
Each point represents one galaxy snapshot. The snapshots from 
the different galaxy simulations are colour-coded as indicated in 
the legend. The best fit to the data is shown as a dotted line, 
while the typical error on each data point is indicated in the 
bottom left corner. Top; showing the measured (unaltered) slopes 
of the cluster age distributions. Bottom: showing slopes that are 
corrected for the variation of the star formation rate (SFR). 

physical meaning because we determine ph.spR for a sphere 
of which a non-negligible fraction is constituted by empty 
space. If the slopes of the age distributions are corrected 
for the variation of the SFR instead of using the raw age 
distributions, we obtain the relation 

/3 = C- 0.68 log ph.SFR, (29) 

with C = —5.04. The errors on the fitted slopes in Eqs. 28 
and 29 are smaller than the listed accuracy. The fitted slopes 
vary by less than 0.03 if the galaxy in the top-left corner of 
both panels in Fig. 11 (IdG) is excluded, which underlines 
the reliability of the fits. 

The physical correlation between the slope of the age 
distribution and the star formation rate density is best de- 
scribed by Eq. 29, because in isolated discs [3 is independent 
of the variation of the SFR. Conversely, the relation between 
a and ph.SFR (Eq. 28) would be relevant for comparison 
with observations. Either way, the implication of both re- 
lations is that the rate of cluster disruption increases with 
the star formation rate density. In Kruijssen et al. (2010), 



5 GALAXY MERGERS 

We now consider the galaxy merger simulations from Ta- 
ble 2. We discuss the evolution of individual clusters, as 
well as the evolution of the cluster population as a whole. 
The section is concluded with a discussion of the cluster 
population in a merger remnant. 

5.1 The evolution of individual clusters in galaxy 
mergers 

Similar to Fig. 5 for disc galaxies, the evolution of the orbits 
and masses of three 'representative' clusters from simulation 
lm2 are shown in Fig. 12. As in Sect. 4.1, the clusters have 
comparable initial masses {Mi ~ 1.5 x 10* M©) and times 
of formation {t ~ 0.12 Gyr), and the differences in their 
evolution are the result of their contrasting orbits in different 
environments. 

The snapshots in Fig. 12 follow the merger during the 
first pericentre passage, when the orbital differences between 
the clusters are partially conserved. This is not the case dur- 
ing the final coalescence of the two galaxies, when violent re- 
laxation randomises the cluster orbits. Just like in isolated 
disc galaxies (Fig. 5), the cluster closest to the centre of the 
galaxy has a low survival chance and is disrupted within 
~ 200 Myr. The two other clusters survive the first passage 

The starbursts in galaxy mergers are characterised by high 
star formation rate densities, yet the slope of the cluster age dis- 
tribution is reported to be —1 (Whitmore et al. 2007), seemingly 
contradicting Fig. 11 and Eq. 28. We discuss the inclusion of 
galaxy mergers in Sect. 5. 
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Figure 12. Evolution of the orbits and masses of tliree clusters in galaxy merger simulation lm2 during the first pericentre passage 
of the galaxies. From top to bottom, the consecutive panels show the situation at different times t, while from left to right the orbital 
evolution in the x-y plane (face-on), the orbital evolution in the x-z plane (edge-on), and the mass evolution are shown. The respective 
clusters are represented by dark red solid lines, red dotted lines, and blue dashed lines. If at any particular snapshot a cluster is still 
undisrupted, its position and mass are marked with thick dots. The orbital trajectories remain visible after the clusters are disrupted. 
The small grey dots in the x-y and x-z plane views map the distribution of the gas particles in the simulation. 
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of the galaxies and experience different evolutionary scenar- 
ios. One is ejected from the disc of its parent galaxy (the red 
dotted cluster in Fig. 12) , together with all the surrounding 
gas and stars, and ends up in the trailing tidal tail of the 
galaxy. It has a low velocity with respect to the tidal tail, 
but it does experience an intermediate tidal shock when en- 
tering the tidal arm at t = 0.32 Gyr, and a strong tidal shock 
when it hits the densest part at f = 0.45 Gyr, leading to the 
disruption of the cluster. The other cluster (blue dashed in 
Fig. 12) is ejected from the disc as well, but it decouples 
from the surrounding gas. This occurs because the gas col- 
lides with the other galaxy and is shocked, which slows it 
down to form the bridge between both galaxies. By contrast, 
the cluster retains a ballistic orbit and becomes part of the 
stellar halo surrounding the galaxies. As a result, the tidal 
field strength decreases and the frequency of tidal shocks be- 
comes low, since the cluster is only shocked twice per orbit. 
The tidal shocks occur when the cluster crosses the bridge 
or the tidal arm and cause it to lose only a few percent of 
its mass. Under these conditions, the expected disruption 
time of the cluster is several gigayears. Even though the 
cluster mass is only 10* M0, this could increase to 10 Gyr 
or more when the tidal arms disperse and the merger con- 
sumes the remaining gas, provided that the cluster does not 
fall back into the central region of the merger. This shows 
that long-lived constituents of the stellar halo surrounding 
giant elliptical galaxies are already produced during the first 
pericentre passage of the progenitor galaxies (see Sect. 5.3). 

The cluster evolution depicted in Fig. 12 illustrates the 
mechanisms of cluster migration and natural selection that 
were explained in Sect. 4.4 and Fig. 10. The cluster that de- 
couples from the gas and is ejected into the stellar halo expe- 
riences a disruption rate that decreases as the cluster ages, 
showing how migration influences the evolution of the clus- 
ter population. On the other hand, the cluster that initially 
resides close to the galactic centre is quickly disrupted by the 
tidal shock of the first pericentre passage, while the two sur- 
viving clusters were situated in less dense environments and 
therefore survive. This shows how natural selection governs 
which clusters survive, and that the mean disruption rate of 
the population decreases with age as clusters in disruptive 
environments are destroyed. 

The mass loss histories of the clusters in Fig. 12 can 
be understood by considering the evolution of the tidal field 
strength and the heating by tidal shocks. Similar to Fig. 6 
in Sect. 4.1, this is shown in Fig. 13 for the clusters in the 
merger. It conflrms that the short-lived cluster indeed ex- 
periences a tidal field strength and tidal shock heating that 
is only rivaled by the cluster that ends up in the halo. The 
reason that the halo cluster is not disrupted on the same 
timescale as the short-lived cluster is that its migration to 
the halo occurs before disruption would have led to its com- 
plete dispersion, thereby decreasing the tidal field strength 
it experiences. The halo cluster therefore only sustains en- 
hanced disruption when it passes through the bridge be- 
tween the two galaxies (at t — 0.5 Gyr), while the short- 
lived cluster stays in a dense environment and is completely 
disrupted by two subsequent tidal shocks. By contrast, the 
cluster in the tidal tail continuously experiences tidal shocks 
and a stronger tidal field than the halo cluster, because it 
is moving with the tidal tail and its environment does not 
change. This leads to an almost constant mass loss rate. 




0.2 0.3 0.4 0.5 
t [Gyr] 

Figure 13. Evolution of the cluster mass and the tidal field for 
the three clusters from Fig. 12, indicated by the same colours 
and line styles. The diamonds in the middle and bottom panel 
indicate the times of each snapshot. Top: The mass evolution. 
Middle: Evolution of the tidal field strength experienced by each 
cluster, defined as the largest eigenvalue of Eq. 8 (see Sect. 2.2.2). 
Bottom: Running integral of the total amount of shock heating 
experienced by the cluster (see Sect. 4.1 and Eq. 25). 



which is enhanced by the tidal shocks occurring when the 
cluster first enters the tidal tail and also when it hits the 
dense centre of the tail. This second tidal shock occurs in 
between two snapshots and the corresponding shock heat- 
ing is therefore not visible in Fig. 13. The decrease of the 
mean tidal field strength and tidal shock heating with age 
illustrate the mechanism of natural selection, i.e. the higher 
survival chances of clusters in quiescent tidal environments. 
The effects of cluster migration and natural selection are 
stronger in galax;y mergers than in isolated disc galaxies, 
because both mechanisms are driven by the variation of the 
environment with time and space. Such variations are evi- 
dently more common in galaxy mergers than in disc galaxies. 



Star cluster populations in galaxy simulations 21 



100 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 
t [Gyr] 

Figure 14. Time evolution of (top) tlie star formation rate and 
(bottom) tlic fitted slope of the age distribution in the range 
log (r/yr) = 7.7-9 for merger simulation lml4, with the red solid 
line denoting the fit to the actual age distribution a, and the blue 
dotted line denoting the fit to the SFR-corrected age distribution 
/3. The dashed vertical lines indicate the moments of first and 
second pericentre passage, and the shaded area marks the period 
over which the final coalescence occurs. 

5.2 The age distributions of star clusters in 
galaxy mergers 

The variation of the environment in galaxy mergers leads 
to a corresponding evolution of the cluster age distribution. 
Similar to Sect. 4.4, we have fitted the slope of the clus- 
ter age distributions in the range log(T/yr) = 7.7-9 for 
all galaxy merger simulations, up to the moment of their 
largest starburst, which typically occurs early on during the 
final coalescence of both galaxies. The slope is not fitted for 
later times, because the gas is rapidly consumed during the 
starburst. At first, this makes the variation of the cluster 
formation rate dominate the shape of the cluster age dis- 
tribution, implying that a power law fit is very inaccurate, 
while later on the age distribution becomes discontinuous 
due to episodes without any surviving clusters (see Sect. 5.3 
and Fig. 17). Similar to Sect. 4.4, we consider all clusters 
when constructing the age distribution, i.e. the samples are 
mass-limited with M ^ 100 Mq . 

In Fig. 14, the star formation history and evolution of 
the fitted slope of the age distribution are shown for merger 
simulation lml4 (see Table 2). The slope widely changes 
over the course of the merger, and is shallowest at the times 
when the star formation rate and star formation rate den- 



sity are highest, with typical slopes between —0.5 and — 1. 
This behaviour is opposite to what is found in Sect. 4.4 for 
isolated disc galaxies, in which the age distribution becomes 
steeper for higher star formation rate densities. As was dis- 
cussed in Sect. 4.4, a shallow slope indicates that cluster mi- 
gration and natural selection are important, i.e. that there 
are large density contrasts in a galaxy, particularly between 
star forming regions and their surroundings. In isolated disc 
galaxies, such a large contrast exists for galaxies with an 
overall low gas density, which then contrasts with the dense 
star forming regions. This low gas density translates to a low 
star formation rate density, and gives the relation of Eqs. 28 
and 29. In galaxy mergers, the effects of cluster migration 
and natural selection are largest at the height of the interac- 
tion. At that point, the star formation rate (density) peaks, 
because all gas is funneled to the central regions, leading to a 
pronounced density contrast between the concentrated star 
forming volume and the surrounding regions, which hardly 
contain any gas. In the meanwhile, the ongoing interaction 
ejects the clusters into the gas-poor stellar halo. The result 
is visible in Fig. 14, in which the slope of the age distribu- 
tion evolves to shallower slopes during the starbursts. The 
extreme slopes in between the starbursts are typically —2.5 
to —3, which is steeper than found in isolated discs. The 
reason is illustrated below, in the discussion of Fig. 15. 

Another interesting feature of Fig. 14 is the difference 
between the actual slope a and the SFR-corrected slope (3. 
Because for /? the variation of the SFR is divided out, one 
would expect it to have a more stable evolution than a. 
However, this is not the case in Fig. 14, where the varia- 
tion of the SFR-corrected slope is larger than that of the 
actual slope. This is the result of the mechanism identified 
in Kruijssen et al. (2010), who find that the gas density in 
starbursts is so high that the young clusters formed in the 
starburst are disrupted on much shorter timescales than in 
isolated galaxies, even to the extent that the total number 
of star clusters decreases during a starburst. This counter- 
intuitive result is mainly due to the lowest mass clusters, 
which are the most numerous for a power law initial mass 
function with a negative slope^^. This large number of low 
mass clusters is susceptible to disruption by the strong tidal 
shocks in a starburst region. The surprising consequence is 
that after a certain time interval, the age distribution of all 
clusters lacks clusters in the age range corresponding to the 
starburst. After the starburst, the peak in the cluster age 
distribution shifts to ages just before the maximum of the 
starburst (also see Sect. 5.3 and Fig. 18), when the clusters 
are still formed in a less violent setting than at the height of 
the burst, and can be ejected from their primordial regions 
before the starburst reaches its maximum (like the halo clus- 
ter of Fig. 12). 

The evolution of the age distribution is compared to 
the star formation history in Fig. 15, which shows the evo- 
lution of the age distribution at several times after the ma- 
jor starburst in simulation lml4. It illustrates several of 
the points from the previous paragraphs. The first age dis- 
tribution (at t = 2.02 Gyr) shows the cause of the steep 



The index —2 of the cluster initial mass function adopted in 
this study implies that every decade in cluster mass initially has 
ten times more clusters than the next decade. 
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Figure 15. Time-evolution of tlie cluster age distribution (solid 
lines) and star formation history (dotted lines) shortly af- 
ter the second passage of merger simulation lml4 (at t fa 
2 Gyr). Prom top to bottom, the distributions are shown at 
times t = {2.02,2.07,2.18,2.41} Gyr, corresponding to about 
{0,50,150,400} Myr after the starburst. For each line, the mo- 
ment of the starburst is marked with a diamond, while the peak 
in the cluster age distribution is indicated with a triangle. Each 
age distribution is shifted down by 5 dex with respect to the dis- 
tribution above it. The star formation histories are normalised to 
match the corresponding age distribution at the left end of the 
lines. For each pair of distributions, the age-offset between the 
peaJjs At is indicated. 



slope of about —3 just before the second starburst. The fit- 
ted slope has steepened relative to isolated galaxies (com- 
pare Fig. 11) due to a deficit of clusters at ages close to 
r = 1 Gyr, which corresponds to the first starburst, when 
the high densities triggered enhanced cluster disruption. The 
same mechanism causes an age-offset between the moment 
of the second starburst and the peak in the age distribu- 
tion, which first emerges when the clusters formed in the 
starburst have had the time to be disrupted by their envi- 
ronment. This disruption is evident from the minimum in 
the age distribution at ages slightly younger than the star- 
burst. The exact moment when the offset between the peaks 
becomes visible depends on the duration and strength of the 
starburst, but it typically appears 100 Myr after the star- 
burst. The offset grows from At = 14.5 Myr at t = 2.07 Gyr 
to At = 132 Myr at t = 2.41 Gyr. As shown in Fig. 15, it is 
best seen about 150 Myr after the burst. When considering 
only the massive clusters (M > 10* Mq), which are much 



less numerous than the low-mass clusters, the deficit of clus- 
ters is less prominent. In the extreme case, the offset of the 
peak in the cluster age distribution with respect to the mo- 
ment of maximum star formation corresponds to the time 
interval between the onset and the peak of the starburst. 

The age-offset between the starburst and the peak in 
the cluster age distribution has an interesting consequence 
in relation to Fig. 14. When dividing the cluster age dis- 
tribution by the star formation history for a galaxy merger 
with a recent starburst, the age range corresponding to the 
starburst will contain even fewer clusters than without the 
correction for the SFR. As a result, the variation of the age 
distribution is enhanced with respect to the actual ago dis- 
tribution. This causes the larger variation of P than that of a 
in Fig. 14. The offset between the peaks in the age distribu- 
tions of the clusters and stars is also seen when considering 
the formation history of the clusters that survive the merger 
(sec Sect. 5.3 and Fig. 18), which shows that these clusters 
arc typically formed before instead of during the starburst 
maximum. It depends on the accuracy of the age determina- 
tions of real clusters whether the offset can be distinguished 
observationally, especially because it is less pronounced for 
the high cluster masses to which observations are naturally 
limited. 

In order to consider the relation between the slope of 

the age distribution and the star formation rate density, we 
have used the same approach as in Sect. 4.4 to determine a 
measure of the star formation rate density in galaxy merg- 
ers. For both galaxies, we determine the half-mass radius of 
the gas distribution and add the enclosed volumes, leaving 
out any overlap between both spheres. To avoid artificially 
low star formation rate densities, the tidal arms are omit- 
ted when calculating the half-mass radius by neglecting all 
material beyond 100 kpc from the centre of mass of the sim- 
ulation. The plane of the fitted age distribution slope versus 
star formation rate density is shown in Fig. 16 for all galaxy 
merger simulations, also including the data from the galaxy 
disc simulations (see Fig. 11). As explained above in the 
discussion of Fig. 14, the galaxy mergers do not follow the 
relation between slope and star formation rate density that 
holds for isolated disc galaxies. Instead, during starbursts 
they typically move to shallower slopes and higher star for- 
mation rate densities, i.e. up and to the right in Fig. 16. 
The large scatter on the points of the galaxy merger sim- 
ulations arises because of the wide range of possible age 
distribution slopes over the course of a single merger, which 
is also present in Fig. 14. The scatter is also increased by our 
method of estimating a measure for the star formation rate 
density, which only allows for an order-of-magnitude anal- 
ysis because it is sensitive to the global dynamical changes 
during the merger. 

The typical evolution of a galaxy merger in the diagram 
of Fig. 16 is illustrated by the evolutionary track of simula- 
tion lml4, which goes through three phases. Initially, both 
galaxies reside on the relation for isolated disc galaxies (dot- 
ted line and cross). For simulation lml4, this is not shown 
in Fig. 16, because it occurs too early on in the simulation 
and insufficient clusters exist in the fitted age range. The 
evolutionary track starts at the top middle of the diagram, 
during the first pericentre passage, when the star formation 
rate density is still intermediate (ph,SFR ~ 10~* Mq kpc~^) 
and cluster migration and natural selection are important. 
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Figure 16. Relation between the fitted logarithmic slope of the 
cluster age distribution in the age range log(r/yr) = 7.7-9 and 
the mean star formation rate density pjj sfr. Each point repre- 
sents one galaxy snapshot. The isolated disc galaxies are shown 
as light grey points, and the galaxy mergers are shown as red 
points. As in Fig. 11, the fit to the isolated disc galaxies is repre- 
sented by a dotted line, while the typical error on each data point 
is indicated in the bottom left corner. Top: the measured (unal- 
tered) slopes of the cluster age distributions. Bottom: slopes that 
are corrected for the variation of the star formation rate (SFR). 
The solid line in both panels indicates the evolutionary track of 
simulation lml4, of which the evolution of the slope is shown in 
Fig. 14. The mean slope and Ph.SFR of the progenitor galaxies 
(IdB) is indicated with a cross. 

resulting in a shallow age distribution. In between both peri- 
centre passages, it returns to the relation for isolated discs 
because the discs evolve back to a quasi-isolated state as in 
Fig. 11, but with a slightly higher star formation rate den- 
sity. This changes just before the final coalescence, when the 
density contrast between the starburst region and the sur- 
roundings becomes important again, moving the galaxy to 
the top right of Fig. 16. 

5.3 The cluster population of merger remnants 

After a galaxy merger is completed and both galaxies have 
transformed into a single elliptical or SO galaxy, the forma- 
tion of stars and clusters ceases or proceeds at a low rate 
(SFR < 0.5 M0 yr~^). As a result, the vast majority of 
clusters in a merger remnant is old, with ages dating back 
to the first and second pericentre passages of the interaction. 
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Figure 17. Logarithmic age distribution dA'^/d log (r/yr) of the 
clusters with ages r ^ 1 Gyr in the merger remnant of simula- 
tion lml4, shown for the snapshot at i = 4.9 Gyr. The vertical 
dashed lines indicate the moments of first (right) and second (left) 
pericentre passage, while the shaded area marks the period over 
which the final coalescence occurs. 

A first indication of when and where these clusters (the 'sur- 
vivors') were formed is obtained from their age distribution, 
which is shown in Fig. 17 for the cluster population older 
than 1 Gyr of simulation lml4. The age distribution shows 
that most of the survivors are formed approximately at the 
times of the first and second pericentre passages, just be- 
fore the corresponding starbursts. During the last part of 
the coalescence, some more clusters are formed that survive 
the merger. Interestingly, no clusters with ages correspond- 
ing to the onset of the coalescence exist in the merger rem- 
nant, because the violent gas influx and the resulting high 
gas density shortens the lifetimes of the clusters that are 
formed under these conditions. 

A more precise picture of the origin of the cluster pop- 
ulation in the merger remnant is obtained by considering 
their cumulative formation history and the radial evolution 
of their population. This is shown in Fig. 18, which follows 
the time evolution of the (cumulative) relative formation his- 
tory and the half-number radius for three groups of objects: 
all survivors formed since the start of the simulation (giving 
a cumulative fraction) , the survivors formed during the last 
200 Myr, and all star particles formed since the start of the 
simulation (also giving a cumulative fraction). Contrary to 
the half-mass radius of the gas in Sects. 4.4 and 5.2, the half- 
number radius considered here is not defined with respect 
to the centre of the appropriate galaxy, but with respect to 
the centre of mass of the entire simulation. 

The cumulative formation history of the survivors shows 
that each of both pericentre passages contributes about 30- 
60% of the old cluster population in the merger remnant. 
The precise distribution of percentages depends on the or- 
bital geometry of the merger and on the properties of the 
progenitor galaxies. In the case of simulation lml4, which is 
shown in Fig. 18, the galaxies pass each other on near-polar 
orbits, yielding a weaker starburst than a head-on or co- 
planar encounter and leaving some gas for post-merger star 
formation. For more violent starbursts, all gas is consumed 
and no young clusters exist in the merger remnant. 
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Figure 18. (Cumulative) formation liistory and radial evolution 
of the clusters that will survive the galaxy merger of simulation 
lml4, i.e. those that are still present at t = 4.9 Gyr. Top: For each 
time t, the figure shows the fraction of the surviving cluster popu- 
lation that has been formed since the start of the simulation (red 
solid line) and the fraction that was formed during the 200 Myr 
preceding t (blue dotted line). The dark red dashed line shows 
the cumulative fraction of star particles that have been formed 
since the start of the simulation. Bottom: Half-number radius of 
all present survivors (red solid line), of the survivors that were 
formed during the 200 Myr interval before time t (blue dotted 
line), and of the star particles that have been formed since the 
start of the simulation (dark red dashed line). Stars and clusters 
formed in the range t = 4—5 Gyr are ignored. 



enhanced star formation rate before the ejection and the 
increased survival chances of halo clusters implies that the 
ejected clusters constitute a large part of the post-merger 
population of survivors. 

The ejection of clusters can also be seen by consider- 
ing the half-number radii of the system of (recently formed) 
surviving clusters and of the stars in Fig. 18. The pericen- 
tre passages of the two galaxies are visible as minima in 
the evolution of the half-number radius of the star particles. 
Already during the first passage, the half-number radius of 
the clusters exceeds that of the star particles, because the 
clusters that are ejected from the discs of both galaxies have 
higher survival chances than the clusters that stay confined 
to the discs. This effect becomes even more important dur- 
ing the second passage and final coalescence of the galaxies, 
during which the half-number radius of the clusters hardly 
changes, but the star particles end up in a much smaller vol- 
ume. While this could suggest that almost no survivors are 
formed at small radii, the half-number radius of the recently 
formed surviving clusters proves the contrary. During and 
shortly (~ 50 Myr) after the second pericentre passage, the 
spatial distribution of the recently formed survivors (blue 
dotted line in Fig. 18) is as confined as the spatial distribu- 
tion of star particles. This illustrates that the clusters may 
be formed in the galactic discs, but are subsequently ejected 
due to the dynamical interaction of the galaxies, increasing 
their chances for survival. At later times (> 50 Myr after 
the pericentre passage), the survivors are formed at differ- 
ent locations than the stars, because the clusters that are 
formed in the centre of the starburst are disrupted. These 
two examples of natural selection imply that the spatial dis- 
tribution of the star cluster population in merger remnants 
does not follow the distribution of the stars, but is spatially 
more extended. 



6 DISCUSSION 

In this section, we provide a summary and a discussion of 
the possible improvements and potential applications of our 
method. 



The assembly history of the stellar mass is distributed 
over both pericentre passages in a way that is similar to that 
of the clusters, even though the first passage gives rise to a 
much smaller starburst than the second passage. The stellar 
mass formed in both passages is comparable because the du- 
ration of the first starburst exceeds that of the second. The 
most remarkable difference between the formation history 
of the star particles and the surviving clusters is that they 
are offset with respect to each other. The surviving clusters 
are typically formed at earlier times than the star particles, 
which was also mentioned in Sect. 5.2 and the discussion of 
Fig. 17. Most of these survivors were ejected into the stel- 
lar halo during the pericentre passages and survived because 
halo clusters experience a lower disruption rate than clusters 
residing in the discs of both galaxies. These ejected clusters 
were formed before the starburst, because the onset of ejec- 
tion into the halo precedes the moment of peak starburst 
intensity by ~ 200 Myr. The combination of an already 



6.1 Summary 

We have presented numerical simulations of isolated and 
merging disc galaxies, in which a sub-grid model for the for- 
mation and evolution of the entire star cluster population 
is included. The description for the star clusters is semi- 
analytic and includes a model for their internal dynamical 
evolution and the resulting changes of the stellar mass func- 
tion within the clusters. The prescription for cluster dis- 
ruption has been validated by comparing to A'^-body sim- 
ulations of dissolving star clusters, giving good agreement. 
When considering individual clusters within our simulations, 
the tidal field strength and tidal shocks are found to have 
a clear effect on the mass loss histories of the clusters. This 
provides a verification of the presented method. 

One of the key advantages of the model is that it shows 
how the disruption rate of clusters varies in time and space. 
We have used our disc galaxy simulations to assess the im- 
plications of this for characteristic properties of the cluster 
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populations. We find that the mean age of the cluster popu- 
lation increases with galactocentric radius, because the dis- 
ruption rate and the cluster formation rate are highest near 
the galactic centre. This is also found in observations of the 
clusters in M51 (Giclos ct al. 2005) and the Milky Way (van 
den Bergh & McClure 1980; Froebrich et al. 2010). The rel- 
ative contribution of tidal shocks to the disruption of star 
clusters is found to be ~ 80%, which weakly increases with 
increasing gas fraction of the galactic disc. A similar value 
was found by Lamcrs & Gieles (2006) from an analysis of 
clusters in the solar neighbourhood. 

The combination of disruption due to two-body relax- 
ation, tidal shocks, and their variation in time and space 
affects the slope of the cluster age distribution through two 
main mechanisms that lead to the same result. 'Cluster mi- 
gration', i.e. the motion of clusters away from their forma- 
tion sites, and 'natural selection', i.e. the higher survival 
probability of clusters in quiescent environments, both im- 
ply that the mean disruption rate decreases with age. In 
the extreme case, this can cause an age distribution with a 
single logarithmic slope of —1 over the majority of the age 
range, instead of the classical flat distribution at young ages 
combined with a steep decline at old ages. For isolated disc 
galaxies, the effects of cluster migration and natural selec- 
tion are largest in low gas density galaxies, because these 
have higher gas density contrasts between star forming re- 
gions and their surroundings. Combining this with the re- 
lation between gas density and star formation rate density 
(Schmidt 1959; Kennicutt 1998b), we obtain a clear correla- 
tion between the star formation rate density and the slope 
of the disruptive (old) end of the age distribution, which is 
steeper for higher star formation rate densities. 

Our simulations of galaxy mergers show that the dis- 
ruption rates of clusters vary widely and depend on their 
orbital histories during the merger. The clusters that re- 
side in the central regions of the galaxies arc disrupted on 
short timescales, while clusters that are ejected into the stel- 
lar halo can survive for several gigayears. The mechanisms 
of cluster migration and natural selection axe prevalent in 
galaxy mergers, because the environment of clusters strongly 
varies in time and space. As a result, the fitted slope of the 
cluster age distribution (in the range log(T/yr) = 7.7-9) 
evolves from —0.5 or —1 during the starbursts, when the 
contrast between the concentrated star forming volume and 
its surroundings is largest, to —2.5 or —3 in between the 
periccntre passages, when the discs evolve back to a quasi- 
isolated state. This is a fundamental physical difference com- 
pared to isolated galaxies, in which the density contrast be- 
tween star forming regions and their surroundings is largest 
for galaxies with low star formation rate densities. 

The star clusters that survive the merger and populate 
the merger remnant are typically formed at the moments 
of the pericentre passages, i.e. slightly before the starbursts 
that occur during a galaxy merger. These clusters consti- 
tute a large fraction (30-60% per pericentre passage) of the 
survivors for two reasons. Firstly, they are formed in large 
numbers, because the star formation rate already increases 
before the peak of the starburst. Secondly, during the peri- 
centre passage, the formed clusters are ejected into the stel- 
lar halo, where the disruption rate is low and the survival 
chance is high. The clusters that are produced in the central 
region during the peak of the starburst are short-lived and 



disrupt before they can migrate to the halo. As a result, a 
peak in the star formation rate docs not necessarily corre- 
spond to a peak in the cluster age distribution. Depending 
on the properties of the starburst and the time that elapsed 
since it occurred, both peaks will be offset with respect to 
each other. 

This paper shows that the variability of the disruption 
rate in time and space has a pronounced impact on the prop- 
erties of cluster populations in a range of galactic environ- 
ments. It affects the spatial distribution of clusters, their age 
distribution, and the evolutionary histories of the clusters 
that survive until the present day. As discussed in Sect. 1, 
it has been common practice in literature to adopt a single, 
"mean" disruption rate for the entire cluster population of 
a galaxy. While this approach holds many advantages due 
to its simplicity, we now see that the resulting cluster pop- 
ulations have very different properties than those ensuing 
from a more realistic setting, in which the effects of the for- 
mation, disruption, and orbital histories of the clusters are 
intertwined. 



6.2 Improvements 

While the presented model gives a more detailed description 
of the formation and evolution of star cluster populations 
than before, there are several points at which it could be 
improved. We discuss five key improvements. 

(1) The current treatment for star formation uses one gas 
particle per spawned star particle. This implies that the par- 
ticle mass limits the maximum cluster mass, because it is not 
possible to form clusters that are more massive than the star 
particle which they are part of (see Sect. 2.2.1). As a result, it 
is not beneficial to increase the resolution of the simulation, 
because it will decrease the maximum cluster mass below 
the current value of ~ lO'"''* M©. Especially when consid- 
ering galaxy mergers, in which clusters with masses around 
10^ Mq should be produced, improving this would be very 
relevant. We intend to include a group-finding algorithm in 
the near future, which will evaluate the Jeans criterion for 
groups of gas particles. This would enable the formation of a 
single star particle out of multiple gas particles, and will also 
allow us to increase the resolution of the simulations with- 
out compromising the mass range of the cluster population. 
In addition to giving a more realistic description of the star 
formation process, this would also enable us to resolve the 
ISM down to smaller scales, and improve the description of 
cluster disruption due to tidal shocks. 

(2) Supermassive black holes (SMBHs) and the possible 
feedback from SMBHs are presently not included. The vast 
majority of star clusters resides in the range where the tidal 
field due to the SMBH can be neglected, so the disruption 
rate of star clusters is not directly affected by the omission of 
SMBHs. An indirect effect of the presence of SMBHs could 
be important in galaxy mergers, during which feedback from 
SMBHs may be responsible for the expulsion of all gas from 
the galaxy (Di Matteo et al. 2005). This would disrupt any 
gas discs that may reform in the merger remnant and would 
halt further formation of star and clusters. Because it is a 
second order effect for the problem we are addressing, and 
because there are currently no definitive models for SMBH 
feedback (Pelupessy 2007; Sijacki et al. 2010), we have cho- 
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sen to omit SMBHs in the present model. Whenever a more 
conclusive model for SMBH feedback becomes available, it 
will be included in our model. 

(3) We have approximated the evolution of the half-mass 
radius of star clusters with a simple power law dependence 
on the cluster mass, fixing the normalisation and power law 
index by means of a comparison to iV-body simulations of 
dissolving clusters on eccentric orbits. This is important, be- 
cause the disruption timescale due to tidal shocks depends 
on the half-mass density. Using the adopted relation, we 
reproduce the disruption times found in the A'^-body sim- 
ulations. Even though the relation is consistent with the 
theoretically expected relation in the 'mass loss-dominated' 
regime from Gieles et al. (2011), a better approach would be 
to adopt a prescription for the half-mass radius that has a 
more extensive physical foundation. Unfortunately, current 
mass-radius relations in literature are based on the evolu- 
tion of clusters in a smooth galactic potential, and depend 
on the galactocentric radius (Gieles et al. 2011). While this 
is accurate for globular clusters on orbits with a low eccen- 
tricity, it does not work for clusters orbiting within a galactic 
disc or in galaxy mergers, where the tidal field is erratic due 
to the non-uniform distribution of the gas. An appropriate 
model for the evolution of the half-mass radius in such an 
environment could be obtained by feeding an erratic tidal 
field into A'^-body simulations of star clusters and monitor- 
ing their structural evolution. Such an analysis is well be- 
yond the scope of the present work, and we will update the 
mass-radius relation whenever a better description becomes 
available. 

(4) At present, the model does not include a description for 
chemical enrichment, and consequently all clusters have the 
same metallicity. While this has a negligible effect on the 
mass evolution of the clusters, their photometry is affected 
(see Kruijssen & Lamers 2008 for a quantitative analysis). 
Moreover, including a prescription for the chemical evolution 
of the star cluster population would enable a better compar- 
ison with (spectroscopic) observations, in which chemical 
abundances can be established with a generally higher ac- 
curacy than other properties such as cluster ages. It would 
also allow us to investigate the relation between metallicity 
and other characteristics of the cluster population, and to 
improve the model for star formation, which depends on the 
chemical composition of the gas. We aim to include a model 
for chemical enrichment in a future work. 

(5) The cluster formation efficiency (CFE), i.e. the frac- 
tion of stars that is formed in a clustered form, is assumed 

to be constant. This implies that the exact value acts as a 
normalisation of the total number of clusters, leaving it as 
a free parameter. It is set to 90% to obtain better statis- 
tics for the simulated cluster populations (see Sect. 2.2.1). 
However, there have been suggestions that the CFE depends 
on the local environment, particularly on quantities like the 
star formation rate density (Goddard et al. 2010). The exact 
dependence of the CFE is still far from certain, but if there 
exists an environmental dependence, this would affect the 
cluster population by favouring the formation of clusters in 
certain parts of a galaxy. This could also have a secondary 
effect on the cluster population, because cluster disruption 
may also proceed differently in parts with an enhanced CFE. 
Again, an environmental dependence of the CFE will be in- 



cluded when it is better constrained, either from models or 
observations. 

Apart from these main areas for improvement, we will 
keep updating the models as A'^-body simulations and ob- 
servations of clusters in a broader range of environments 
become available. 

6.3 Applications 

In order to trace the formation and evolution of galaxies 
using star cluster populations, it is necessary to investigate 
how different galactic environments affect the cluster popu- 
lation. Our model is a very suitable tool to gain more insight 
into this question, because it relates the evolution of each 
cluster to its (time-dependent) local environment. This im- 
plies a certain flexibility that allows us to apply the model 
to a broad range of galaxies. While a first analysis of the in- 
terplay between galaxies and their star cluster populations 
is already given in this paper, there are many more observ- 
ables of the cluster population that should be investigated 
under different galactic conditions. 

It would be particularly useful to understand the im- 
pact of galaxy mergers on cluster populations, because such 
an understanding enables the use of cluster populations 
to probe merger histories and the hierarchical assembly of 
galaxies. Mergers are recognized as important drivers of star- 
bursts and corresponding cluster formation, which arc fu- 
elled by high gas densities. However, as is shown in this pa- 
per, a high gas density also implies a large disruption rate. It 
is not trivial to determine whether formation or destruction 
dominates. We have considered this question in Kruijssen 
et al. (2010) as a first application of the model, and find that 
the total number of clusters decreases during a merger, be- 
cause the large gas densities result in more destruction than 
formation. Destruction is most prominent for the numerous 
clusters with low masses, whereas for the fewer massive clus- 
ters formation does dominate during certain episodes of the 
galax;y interaction. The corresponding change of the cluster 
mass function could be used as a tracer of the merger type. 

By modeling specific, real galaxies, it is possible to ex- 
plain observed properties of the cluster population and to 
predict its features that presently fall below the detection 
limit. Such case studies will also verify the model, and pos- 
sibly provide constraints on aspects of the model that are 
currently uncertain (sec Sect. 6.2). For instance, by compar- 
ing the observed and modeled star formation rates and the 
number of clusters within a certain mass and age range, it 
will be possible to infer the cluster formation efficiency in a 
particular galaxy^*. 

As is indicated in Sect. 1, the disruption rate of star 
clusters is commonly assumed to be constant when deriving 
the star formation history (SFH) of a galaxy from its star 
cluster population. Although this approximation is conve- 
nient, the thus obtained SFH will differ from the actual one. 

Before being able to combine different observational studies 
to look at trends of the properties of the cluster population with 
galactic environment, the current dichotomy in literature between 
groups drawing different conclusions from the exact same obser- 
vational data (e.g. Chandar et al. 2006; Gieles et al. 2007b) should 
be settled. No model can bridge such differences. 
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The impact of the disruption time on the inferred SFH was 
recently illustrated by Maschberger & Kroupa (2010, Fig. 4), 
who show that it depends on the adopted disruption rate to 
what extent the gap in the age distribution of clusters in the 
Large Magellanic Cloud is reflected in the inferred SFH. For 
their choice of disruption rates, the SFR in the age range 
corresponding to the age gap varies by about 1.5 dex, re- 
sulting in cases in which the SFR does and does not contain 
the ago gap of the cluster age distribution. Because the con- 
ditions within an evolving galaxy vary widely, the impact 
of the time- and space-variation of the disruption rate are 
likely of the same order of magnitude. It is therefore essen- 
tial to resolve how this variation may affect SFHs that are 
inferred from the star cluster population. 

The formation and evolution of star cluster populations 
are the result of several mechanisms that act simultane- 
ously, such as starbursts, feedback, tidal shocks, two-body 
relaxation, cluster migration, natural selection, and many 
other processes. While certain parts may still be uncertain, 
the current understanding of these mechanisms enables the 
modeling of the cluster population in a way that reflects 
the variability and complex nature of real galactic environ- 
ments. Future applications of the model should therefore 
provide new clues to the (co-)evolution of galaxies and their 
star cluster populations. 
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